diff --git a/README.md b/README.md index 1a5e768..e684d21 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,9 @@ # CSiBORG tools ## :scroll: Short-term TODO -- [x] Add the X-ray clusters -- [x] Match the X-ray clusters to Planck -- [ ] Calculate catalogues for all realisations. +- [ ] Find the distribution of particles in the first snapshot - [ ] Implement Max's halo matching +- [ ] Implement a custom model for matchin galaxies to halos. ## :hourglass: Long-term TODO diff --git a/csiborgtools/__init__.py b/csiborgtools/__init__.py index 94a7917..ed976df 100644 --- a/csiborgtools/__init__.py +++ b/csiborgtools/__init__.py @@ -13,4 +13,4 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -from csiborgtools import (io, match, utils, units, fits) # noqa +from csiborgtools import (read, match, utils, units, fits) # noqa diff --git a/csiborgtools/fits/halo.py b/csiborgtools/fits/halo.py index 296b18c..1f50219 100644 --- a/csiborgtools/fits/halo.py +++ b/csiborgtools/fits/halo.py @@ -22,7 +22,7 @@ from os import remove from warnings import warn from os.path import join from tqdm import trange -from ..io import nparts_to_start_ind +from ..read import nparts_to_start_ind def clump_with_particles(particle_clumps, clumps): diff --git a/csiborgtools/io/readobs.py b/csiborgtools/io/readobs.py deleted file mode 100644 index 7b378aa..0000000 --- a/csiborgtools/io/readobs.py +++ /dev/null @@ -1,263 +0,0 @@ -# Copyright (C) 2022 Richard Stiskalek -# This program is free software; you can redistribute it and/or modify it -# under the terms of the GNU General Public License as published by the -# Free Software Foundation; either version 3 of the License, or (at your -# option) any later version. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General -# Public License for more details. -# -# You should have received a copy of the GNU General Public License along -# with this program; if not, write to the Free Software Foundation, Inc., -# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -""" -Scripts to read in observation. -""" - -import numpy -from astropy.io import fits -from astropy.coordinates import SkyCoord -from astropy import units as u -from ..utils import (add_columns, cols_to_structured) - -F64 = numpy.float64 - - -def read_planck2015(fpath, cosmo, max_comdist=None): - r""" - Read the Planck 2nd Sunyaev-Zeldovich source catalogue [1]. The following - is performed: - - removes clusters without a redshift estimate, - - calculates the comoving distance with the provided cosmology. - - Converts `MSZ` from units of :math:`1e14 M_\odot` to :math:`M_\odot` - - Parameters - ---------- - fpath : str - Path to the source catalogue. - cosmo : `astropy.cosmology` object - The cosmology to calculate cluster comoving distance from redshift and - convert their mass. - max_comdist : float, optional - Maximum comoving distance threshold in units of :math:`\mathrm{Mpc}`. - By default `None` and no threshold is applied. - - Returns - ------- - out : structured array - The catalogue structured array. - - References - ---------- - [1] https://heasarc.gsfc.nasa.gov/W3Browse/all/plancksz2.html - """ - data = fits.open(fpath)[1].data - hdata = 0.7 - # Convert FITS to a structured array - out = numpy.full(data.size, numpy.nan, dtype=data.dtype.descr) - for name in out.dtype.names: - out[name] = data[name] - # Take only clusters with redshifts - out = out[out["REDSHIFT"] >= 0] - # Add comoving distance - dist = cosmo.comoving_distance(out["REDSHIFT"]).value - out = add_columns(out, dist, "COMDIST") - # Convert masses - for par in ("MSZ", "MSZ_ERR_UP", "MSZ_ERR_LOW"): - out[par] *= 1e14 - out[par] *= (hdata / cosmo.h)**2 - # Distance threshold - if max_comdist is not None: - out = out[out["COMDIST"] < max_comdist] - - return out - - -def read_mcxc(fpath, cosmo, max_comdist=None): - r""" - Read the MCXC Meta-Catalog of X-Ray Detected Clusters of Galaxies - catalogue [1], with data description at [2] and download at [3]. - - Note - ---- - The exact mass conversion has non-trivial dependence on :math:`H(z)`, see - [1] for more details. However, this should be negligible. - - Parameters - ---------- - fpath : str - Path to the source catalogue obtained from [3]. Expected to be the fits - file. - cosmo : `astropy.cosmology` object - The cosmology to calculate cluster comoving distance from redshift and - convert their mass. - max_comdist : float, optional - Maximum comoving distance threshold in units of :math:`\mathrm{Mpc}`. - By default `None` and no threshold is applied. - - Returns - ------- - out : structured array - The catalogue structured array. - - References - ---------- - [1] The MCXC: a meta-catalogue of x-ray detected clusters of galaxies - (2011); Piffaretti, R. ; Arnaud, M. ; Pratt, G. W. ; Pointecouteau, - E. ; Melin, J. -B. - [2] https://heasarc.gsfc.nasa.gov/W3Browse/rosat/mcxc.html - [3] https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/534/A109#/article - """ - data = fits.open(fpath)[1].data - hdata = 0.7 # Little h of the catalogue - - cols = [("RAdeg", F64), ("DEdeg", F64), ("z", F64), - ("L500", F64), ("M500", F64), ("R500", F64)] - out = cols_to_structured(data.size, cols) - for col in cols: - par = col[0] - out[par] = data[par] - # Get little h units to match the cosmology - out["L500"] *= (hdata / cosmo.h)**2 - out["M500"] *= (hdata / cosmo.h)**2 - # Get the 10s back in - out["L500"] *= 1e44 # ergs/s - out["M500"] *= 1e14 # Msun - - dist = cosmo.comoving_distance(data["z"]).value - out = add_columns(out, dist, "COMDIST") - out = add_columns(out, data["MCXC"], "name") - - if max_comdist is not None: - out = out[out["COMDIST"] < max_comdist] - - return out - - -def read_2mpp(fpath, dist_cosmo): - """ - Read in the 2M++ galaxy redshift catalogue [1], with the catalogue at [2]. - Removes fake galaxies used to fill the zone of avoidance. Note that in - principle additional care should be taken for calculating the distance - to objects [3]. Currently calculated from the CMB redshift, so some - distance estimates may be negative.. - - Parameters - ---------- - fpath : str - File path to the catalogue. - cosmo : `astropy.cosmology` object - The cosmology to calculate distance from redshift. - - Returns - ------- - out : structured array - The catalogue. - - References - ---------- - [1] The 2M++ galaxy redshift catalogue; Lavaux, Guilhem, Hudson, Michael J. - [2] https://cdsarc.cds.unistra.fr/viz-bin/cat/J/MNRAS/416/2840#/article - [3] Improving NASA/IPAC Extragalactic Database Redshift Calculations - (2021); Anthony Carr and Tamara Davis - """ - from scipy.constants import c - # Read the catalogue and select non-fake galaxies - cat = numpy.genfromtxt(fpath, delimiter="|", ) - cat = cat[cat[:, 12] == 0, :] - - cols = [("RA", F64), ("DEC", F64), ("Ksmag", F64), ("ZCMB", F64), - ("DIST", F64)] - out = cols_to_structured(cat.shape[0], cols) - out["RA"] = cat[:, 1] - out["DEC"] = cat[:, 2] - out["Ksmag"] = cat[:, 5] - out["ZCMB"] = cat[:, 7] / (c * 1e-3) - out["DIST"] = cat[:, 7] / dist_cosmo.H0 - return out - - -def read_2mpp_groups(fpath, dist_cosmo): - """ - Read in the 2M++ galaxy group catalogue [1], with the catalogue at [2]. - Note that the same caveats apply to the distance calculation as in - py:function:`read_2mpp`. See that function for more details. - - Parameters - ---------- - fpath : str - File path to the catalogue. - cosmo : `astropy.cosmology` object - The cosmology to calculate distance from redshift. - - Returns - ------- - out : structured array - The catalogue. - - References - ---------- - [1] The 2M++ galaxy redshift catalogue; Lavaux, Guilhem, Hudson, Michael J. - [2] https://cdsarc.cds.unistra.fr/viz-bin/cat/J/MNRAS/416/2840#/article - [3] Improving NASA/IPAC Extragalactic Database Redshift Calculations - (2021); Anthony Carr and Tamara Davis - """ - cat = numpy.genfromtxt(fpath, delimiter="|", ) - - cols = [("RA", F64), ("DEC", F64), ("K2mag", F64), ("Rich", numpy.int64), - ("dist", F64), ("sigma", F64)] - out = cols_to_structured(cat.shape[0], cols) - - out["K2mag"] = cat[:, 3] - out["Rich"] = cat[:, 4] - out["sigma"] = cat[:, 7] - out["dist"] = cat[:, 6] / dist_cosmo.H0 - - # Convert galactic coordinates to RA, dec - glon = cat[:, 1] - glat = cat[:, 2] - coords = SkyCoord(l=glon*u.degree, b=glat*u.degree, frame='galactic') - coords = coords.transform_to("icrs") - out["RA"] = coords.ra - out["DEC"] = coords.dec - - return out - - -def match_planck_to_mcxc(planck, mcxc): - """ - Return the MCXC catalogue indices of the Planck catalogue detections. Finds - the index of the quoted Planck MCXC counterpart in the MCXC array. If not - found throws an error. For this reason it may be better to make sure the - MCXC catalogue reaches further. - - - Parameters - ---------- - planck : structured array - The Planck cluster array. - mcxc : structured array - The MCXC cluster array. - - Returns - ------- - indxs : 1-dimensional array - The array of MCXC indices to match the Planck array. If no counterpart - is found returns `numpy.nan`. - """ - # Planck MCXC need to be decoded to str - planck_names = [name.decode() for name in planck["MCXC"]] - mcxc_names = [name for name in mcxc["name"]] - - indxs = [numpy.nan] * len(planck_names) - for i, name in enumerate(planck_names): - if name == "": - continue - if name in mcxc_names: - indxs[i] = mcxc_names.index(name) - else: - raise ValueError("Planck MCXC identifies `{}` not found in the " - "MCXC catalogue.".format(name)) - return indxs diff --git a/csiborgtools/match/__init__.py b/csiborgtools/match/__init__.py index e09466d..d715465 100644 --- a/csiborgtools/match/__init__.py +++ b/csiborgtools/match/__init__.py @@ -13,6 +13,6 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -from .match import brute_spatial_separation # noqa +from .match import (brute_spatial_separation, RealisationsMatcher) # noqa from .num_density import (binned_counts, number_density) # noqa # from .correlation import (get_randoms_sphere, sphere_angular_tpcf) # noqa diff --git a/csiborgtools/match/match.py b/csiborgtools/match/match.py index c3a01af..ccf8a28 100644 --- a/csiborgtools/match/match.py +++ b/csiborgtools/match/match.py @@ -16,6 +16,7 @@ import numpy from tqdm import tqdm from astropy.coordinates import SkyCoord +from ..read import CombinedHaloCatalogue def brute_spatial_separation(c1, c2, angular=False, N=None, verbose=False): @@ -65,3 +66,130 @@ def brute_spatial_separation(c1, c2, angular=False, N=None, verbose=False): sep[i, :] = dist[sort] return sep, indxs + + +class RealisationsMatcher: + """ + A tool to match halos between IC realisations. Looks for halos 3D space + neighbours in all remaining IC realisations that are within some mass + range of it. + + + Parameters + ---------- + cats : :py:class`csiborgtools.read.CombinedHaloCatalogue` + Combined halo catalogue to search. + """ + _cats = None + + def __init__(self, cats): + self.cats = cats + + @property + def cats(self): + """ + Combined catalogues. + + Returns + ------- + cats : :py:class`csiborgtools.read.CombinedHaloCatalogue` + Combined halo catalogue. + """ + return self._cats + + @cats.setter + def cats(self, cats): + """ + Sets `cats`, ensures that this is a `CombinedHaloCatalogue` object. + """ + if not isinstance(cats, CombinedHaloCatalogue): + raise TypeError("`cats` must be of type `CombinedHaloCatalogue`.") + self._cats = cats + + def search_sim_indices(self, n_sim): + """ + Return indices of all other IC realisations but of `n_sim`. + + Parameters + ---------- + n_sim : int + IC realisation index of `self.cats` to be skipped. + + Returns + ------- + indxs : list of int + The remaining IC indices. + """ + return [i for i in range(self.cats.N) if i != n_sim] + + def cross_knn_position_single(self, n_sim, nmult=5, dlogmass=2): + r""" + Find all neighbours within :math:`n_{\rm mult} R_{200c}` of halos in + the `nsim`th simulation. Also enforces that the neighbours' + :math:`\log M_{200c}` be within `dlogmass` dex. + + Parameters + ---------- + n_sim : int + Index of an IC realisation in `self.cats` whose halos' neighbours + in the remaining simulations to search for. + nmult : float or int, optional + Multiple of :math:`R_{200c}` within which to return neighbours. By + default 5. + dlogmass : float, optional + Tolerance on mass logarithmic mass difference. By default 2 dex. + + Returns + ------- + matches : composite array + Array, indices are `(n_sims - 1, 2, n_halos, n_matches)`. The + 2nd axis is `index` of the neighbouring halo in its catalogue and + `dist`, which is the 3D distance to the halo whose neighbours are + searched. + """ + # R200c, M200c and positions of halos in `n_sim` IC realisation + r200 = self.cats[n_sim]["r200"] + logm200 = numpy.log10(self.cats[n_sim]["m200"]) + pos = self.cats[n_sim].positions + + matches = [None] * (self.cats.N - 1) + # Search for neighbours in the other simulations + for count, i in enumerate(self.search_sim_indices(n_sim)): + dist, indxs = self.cats[i].radius_neigbours(pos, r200 * nmult) + # Get rid of neighbors whose mass is too off + for j, indx in enumerate(indxs): + match_logm200 = numpy.log10(self.cats[i][indx]["m200"]) + mask = numpy.abs(match_logm200 - logm200[j]) < dlogmass + dist[j] = dist[j][mask] + indxs[j] = indx[mask] + # Append as a composite array + matches[count] = numpy.asarray([indxs, dist], dtype=object) + + return numpy.asarray(matches, dtype=object) + + def cross_knn_position_all(self, nmult=5, dlogmass=2): + r""" + Find all neighbours within :math:`n_{\rm mult} R_{200c}` of halos in + all simulations listed in `self.cats`. Also enforces that the + neighbours' :math:`\log M_{200c}` be within `dlogmass` dex. + + Parameters + ---------- + nmult : float or int, optional + Multiple of :math:`R_{200c}` within which to return neighbours. By + default 5. + dlogmass : float, optional + Tolerance on mass logarithmic mass difference. By default 2 dex. + + Returns + ------- + matches : list of composite arrays + List whose length is `n_sims`. For description of its elements see + `self.cross_knn_position_single(...)`. + """ + N = self.cats.N # Number of catalogues + matches = [None] * N + # Loop over each catalogue + for i in range(N): + matches[i] = self.cross_knn_position_single(i, nmult, dlogmass) + return matches diff --git a/csiborgtools/io/__init__.py b/csiborgtools/read/__init__.py similarity index 83% rename from csiborgtools/io/__init__.py rename to csiborgtools/read/__init__.py index 3f22163..d686d5d 100644 --- a/csiborgtools/io/__init__.py +++ b/csiborgtools/read/__init__.py @@ -17,7 +17,7 @@ from .readsim import (get_csiborg_ids, get_sim_path, get_snapshots, # noqa get_snapshot_path, get_maximum_snapshot, read_info, nparts_to_start_ind, # noqa open_particle, open_unbinding, read_particle, # noqa drop_zero_indx, # noqa - read_clumpid, read_clumps, read_mmain, # noqa - merge_mmain_to_clumps) # noqa -from .readobs import (read_planck2015, read_mcxc, read_2mpp, read_2mpp_groups, match_planck_to_mcxc) # noqa + read_clumpid, read_clumps, read_mmain) # noqa +from .make_cat import (HaloCatalogue, CombinedHaloCatalogue) # noqa +from .readobs import (PlanckClusters, MCXCClusters, TwoMPPGalaxies, TwoMPPGroups) # noqa from .outsim import (dump_split, combine_splits) # noqa diff --git a/csiborgtools/read/make_cat.py b/csiborgtools/read/make_cat.py new file mode 100644 index 0000000..6fa5e82 --- /dev/null +++ b/csiborgtools/read/make_cat.py @@ -0,0 +1,331 @@ +# Copyright (C) 2022 Richard Stiskalek, Harry Desmond +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +""" +Functions to read in the particle and clump files. +""" + +import numpy +from os.path import join +from tqdm import trange +from sklearn.neighbors import NearestNeighbors +from .readsim import (get_sim_path, read_mmain, get_csiborg_ids, + get_maximum_snapshot) +from ..utils import (flip_cols, add_columns) +from ..units import (BoxUnits, cartesian_to_radec) + + +class HaloCatalogue: + r""" + Processed halo catalogue, the data should be calculated in `run_fit_halos`. + + Parameters + ---------- + n_sim: int + Initial condition index. + n_snap: int + Snapshot index. + minimum_m500 : float, optional + The minimum :math:`M_{rm 500c} / M_\odot` mass. By default no + threshold. + dumpdir : str, optional + Path to where files from `run_fit_halos` are stored. By default + `/mnt/extraspace/rstiskalek/csiborg/`. + mmain_path : str, optional + Path to where mmain files are stored. By default + `/mnt/zfsusers/hdesmond/Mmain`. + """ + _box = None + _n_sim = None + _n_snap = None + _data = None + _knn = None + _positions = None + + def __init__(self, n_sim, n_snap, minimum_m500=None, + dumpdir="/mnt/extraspace/rstiskalek/csiborg/", + mmain_path="/mnt/zfsusers/hdesmond/Mmain"): + self._box = BoxUnits(n_snap, get_sim_path(n_sim)) + minimum_m500 = 0 if minimum_m500 is None else minimum_m500 + self._set_data(n_sim, n_snap, dumpdir, mmain_path, minimum_m500) + self._nsim = n_sim + self._nsnap = n_snap + # Initialise the KNN + knn = NearestNeighbors() + knn.fit(self.positions) + self._knn = knn + + @property + def data(self): + """ + Halo catalogue. + + Returns + ------- + cat : structured array + Catalogue. + """ + if self._data is None: + raise ValueError("`data` is not set!") + return self._data + + @property + def box(self): + """ + Box object, useful for change of units. + + Returns + ------- + box : :py:class:`csiborgtools.units.BoxUnits` + The box object. + """ + return self._box + + @property + def cosmo(self): + """ + The box cosmology. + + Returns + ------- + cosmo : `astropy` cosmology object + Box cosmology. + """ + return self.box.cosmo + + @property + def n_snap(self): + """ + The snapshot ID. + + Returns + ------- + n_snap : int + Snapshot ID. + """ + return self._n_snap + + @property + def n_sim(self): + """ + The initiali condition (IC) realisation ID. + + Returns + ------- + n_sim : int + The IC ID. + """ + return self._n_sim + + def _set_data(self, n_sim, n_snap, dumpdir, mmain_path, minimum_m500): + """ + Loads the data, merges with mmain, does various coordinate transforms. + """ + # Load the processed data + fname = "ramses_out_{}_{}.npy".format( + str(n_sim).zfill(5), str(n_snap).zfill(5)) + data = numpy.load(join(dumpdir, fname)) + + # Load the mmain file and add it to the data + mmain = read_mmain(n_sim, mmain_path) + data = self.merge_mmain_to_clumps(data, mmain) + flip_cols(data, "peak_x", "peak_z") + + # Cut on number of particles and finite m200 + data = data[(data["npart"] > 100) & numpy.isfinite(data["m200"])] + + # Unit conversion + convert_cols = ["m200", "m500", "totpartmass", "mass_mmain", + "r200", "r500", "Rs", "rho0", + "peak_x", "peak_y", "peak_z"] + data = self.box.convert_from_boxunits(data, convert_cols) + + # Cut on mass. Note that this is in Msun + data = data[data["m500"] > minimum_m500] + + # Now calculate spherical coordinates + d, ra, dec = cartesian_to_radec(data) + data = add_columns(data, [d, ra, dec], ["dist", "ra", "dec"]) + + # Pre-allocate the positions array + self._positions = numpy.vstack( + [data["peak_{}".format(p)] for p in ("x", "y", "z")]).T + + self._data = data + + def merge_mmain_to_clumps(self, clumps, mmain): + """ + Merge columns from the `mmain` files to the `clump` file, matches them + by their halo index while assuming that the indices `index` in both + arrays are sorted. + + Parameters + ---------- + clumps : structured array + Clumps structured array. + mmain : structured array + Parent halo array whose information is to be merged into `clumps`. + + Returns + ------- + out : structured array + Array with added columns. + """ + X = numpy.full((clumps.size, 2), numpy.nan) + # Mask of which clumps have a mmain index + mask = numpy.isin(clumps["index"], mmain["index"]) + + X[mask, 0] = mmain["mass_cl"] + X[mask, 1] = mmain["sub_frac"] + return add_columns(clumps, X, ["mass_mmain", "sub_frac"]) + + @property + def positions(self): + """ + 3D positions of halos. + + Returns + ------- + X : 2-dimensional array + Array of shape `(n_halos, 3)`, where the latter axis represents + `x`, `y` and `z`. + """ + return self._positions + + def radius_neigbours(self, X, radius): + """ + Return sorted nearest neigbours within `radius` or `X`. + + Parameters + ---------- + X : 2-dimensional array + Array of shape `(n_queries, 3)`, where the latter axis represents + `x`, `y` and `z`. + radius : float + Limiting distance of neighbours. + + Returns + ------- + dist : list of 1-dimensional arrays + List of length `n_queries` whose elements are arrays of distances + to the nearest neighbours. + knns : list of 1-dimensional arrays + List of length `n_queries` whose elements are arrays of indices of + nearest neighbours in this catalogue. + """ + if not (X.ndim == 2 and X.shape[1] == 3): + raise TypeError("`X` must be an array of shape `(n_samples, 3)`.") + # Query the KNN + return self._knn.radius_neighbors(X, radius, sort_results=True) + + @property + def keys(self): + """Catalogue keys.""" + return self.data.dtype.names + + def __getitem__(self, key): + return self._data[key] + + +class CombinedHaloCatalogue: + r""" + A combined halo catalogue, containing `HaloCatalogue` for each IC + realisation at the latest redshift. + + Parameters + ---------- + minimum_m500 : float, optional + The minimum :math:`M_{rm 500c} / M_\odot` mass. By default no + threshold. + dumpdir : str, optional + Path to where files from `run_fit_halos` are stored. By default + `/mnt/extraspace/rstiskalek/csiborg/`. + mmain_path : str, optional + Path to where mmain files are stored. By default + `/mnt/zfsusers/hdesmond/Mmain`. + verbose : bool, optional + Verbosity flag for reading the catalogues. + """ + _n_sims = None + _n_snaps = None + _cats = None + + def __init__(self, minimum_m500=None, + dumpdir="/mnt/extraspace/rstiskalek/csiborg/", + mmain_path="/mnt/zfsusers/hdesmond/Mmain", verbose=True): + # Read simulations and their maximum snapshots + # NOTE remove this later and take all cats + self._n_sims = get_csiborg_ids("/mnt/extraspace/hdesmond")[:10] + n_snaps = [get_maximum_snapshot(get_sim_path(i)) for i in self._n_sims] + self._n_snaps = numpy.asanyarray(n_snaps) + + cats = [None] * self.N + for i in trange(self.N) if verbose else range(self.N): + cats[i] = HaloCatalogue(self._n_sims[i], self._n_snaps[i], + minimum_m500, dumpdir, mmain_path) + self._cats = cats + + @property + def N(self): + """ + Number of IC realisations in this combined catalogue. + + Returns + ------- + N : int + Number of catalogues. + """ + return len(self.n_sims) + + @property + def n_sims(self): + """ + IC realisations CSiBORG identifiers. + + Returns + ------- + ids : 1-dimensional array + Array of IDs. + """ + return self._n_sims + + @property + def n_snaps(self): + """ + Snapshot numbers corresponding to `self.n_sims`. + + Returns + ------- + n_snaps : 1-dimensional array + Array of snapshot numbers. + """ + return self._n_snaps + + @property + def cats(self): + """ + Catalogues associated with this object. + + Returns + ------- + cats : list of `HaloCatalogue` + Catalogues. + """ + return self._cats + + def __getitem__(self, n): + if n > self.N: + raise ValueError("Catalogue count is {}, requested catalogue {}." + .format(self.N, n)) + return self.cats[n] diff --git a/csiborgtools/io/outsim.py b/csiborgtools/read/outsim.py similarity index 100% rename from csiborgtools/io/outsim.py rename to csiborgtools/read/outsim.py diff --git a/csiborgtools/read/readobs.py b/csiborgtools/read/readobs.py new file mode 100644 index 0000000..1cfd73c --- /dev/null +++ b/csiborgtools/read/readobs.py @@ -0,0 +1,300 @@ +# Copyright (C) 2022 Richard Stiskalek +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +""" +Scripts to read in observation. +""" + +import numpy +from astropy.io import fits +from astropy.coordinates import SkyCoord +from astropy.cosmology import FlatLambdaCDM +from astropy import units as u +from ..utils import (add_columns, cols_to_structured) + +F64 = numpy.float64 + + +class BaseSurvey: + """ + Base survey class with some methods that are common to all survey classes. + """ + _data = None + _cosmo = None + + @property + def data(self): + """ + Cluster catalogue. + + Returns + ------- + cat : structured array + Catalogue. + """ + if self._data is None: + raise ValueError("`data` is not set!") + return self._data + + @property + def cosmo(self): + """Desired cosmology.""" + if self._cosmo is None: + raise ValueError("`cosmo` is not set!") + return self._cosmo + + @property + def keys(self): + """Catalogue keys.""" + return self.data.dtype.names + + def __getitem__(self, key): + return self._data[key] + + +class PlanckClusters(BaseSurvey): + r""" + Planck 2nd Sunyaev-Zeldovich source catalogue [1]. Automatically removes + clusters without a redshift estimate. + + Parameters + ---------- + fpath : str + Path to the source catalogue. + cosmo : `astropy.cosmology` object, optional + Cosmology to convert masses (particularly :math:`H_0`). By default + `FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728)`. + max_redshift: float, optional + Maximum cluster redshift. By default `None` and no selection is + performed. + + References + ---------- + [1] https://heasarc.gsfc.nasa.gov/W3Browse/all/plancksz2.html + """ + _hdata = 0.7 # little h value of the data + + def __init__(self, fpath, cosmo=None, max_redshift=None): + if cosmo is None: + self._cosmo = FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728) + else: + self._cosmo = cosmo + self.set_data(fpath, max_redshift) + + def set_data(self, fpath, max_redshift=None): + """ + Set the catalogue, loads it and applies a maximum redshift cut. + """ + cat = fits.open(fpath)[1].data + # Convert FITS to a structured array + data = numpy.full(cat.size, numpy.nan, dtype=cat.dtype.descr) + for name in cat.dtype.names: + data[name] = cat[name] + # Take only clusters with redshifts + data = data[data["REDSHIFT"] >= 0] + # Convert masses + for par in ("MSZ", "MSZ_ERR_UP", "MSZ_ERR_LOW"): + data[par] *= 1e14 + data[par] *= (self._hdata / self.cosmo.h)**2 + # Redshift cut + if max_redshift is not None: + data = data["REDSHIFT" <= max_redshift] + self._data = data + + def match_to_mcxc(self, mcxc): + """ + Return the MCXC catalogue indices of the Planck catalogue detections. + Finds the index of the quoted Planck MCXC counterpart in the MCXC + array. If not found throws an error. For this reason it may be better + to make sure the MCXC catalogue reaches further. + + Parameters + ---------- + mcxc : :py:class`MCXCClusters` + MCXC cluster object. + + Returns + ------- + indxs : list of int + Array of MCXC indices to match the Planck array. If no counterpart + is found returns `numpy.nan`. + """ + if not isinstance(mcxc, MCXCClusters): + raise TypeError("`mcxc` must be `MCXCClusters` type.") + + # Planck MCXC need to be decoded to str + planck_names = [name.decode() for name in self["MCXC"]] + mcxc_names = [name for name in mcxc["name"]] + + indxs = [numpy.nan] * len(planck_names) + for i, name in enumerate(planck_names): + if name == "": + continue + if name in mcxc_names: + indxs[i] = mcxc_names.index(name) + else: + raise ValueError("Planck MCXC identifier `{}` not found in " + "the MCXC catalogue.".format(name)) + return indxs + + +class MCXCClusters(BaseSurvey): + r""" + MCXC Meta-Catalog of X-Ray Detected Clusters of Galaxies catalogue [1], + with data description at [2] and download at [3]. + + Note + ---- + The exact mass conversion has non-trivial dependence on :math:`H(z)`, see + [1] for more details. However, this should be negligible. + + Parameters + ---------- + fpath : str + Path to the source catalogue obtained from [3]. Expected to be the fits + file. + cosmo : `astropy.cosmology` object, optional + The cosmology to to convert cluster masses (to first order). By default + `FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728)`. + max_redshift: float, optional + Maximum cluster redshift. By default `None` and no selection is + performed. + + References + ---------- + [1] The MCXC: a meta-catalogue of x-ray detected clusters of galaxies + (2011); Piffaretti, R. ; Arnaud, M. ; Pratt, G. W. ; Pointecouteau, + E. ; Melin, J. -B. + [2] https://heasarc.gsfc.nasa.gov/W3Browse/rosat/mcxc.html + [3] https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/534/A109#/article + """ + _hdata = 0.7 # Little h of the catalogue + + def __init__(self, fpath, cosmo=None, max_redshift=None): + if cosmo is None: + self._cosmo = FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728) + else: + self._cosmo = cosmo + self._set_data(fpath, max_redshift) + + def _set_data(self, fpath, max_redshift): + """ + Set the catalogue, loads it and applies a maximum redshift cut. + """ + cat = fits.open(fpath)[1].data + # Pre-allocate array and extract selected variables + cols = [("RAdeg", F64), ("DEdeg", F64), ("z", F64), + ("L500", F64), ("M500", F64), ("R500", F64)] + data = cols_to_structured(cat.size, cols) + for col in cols: + par = col[0] + data[par] = cat[par] + # Add the cluster names + data = add_columns(data, cat["MCXC"], "name") + + # Get little h units to match the cosmology + data["L500"] *= (self._hdata / self.cosmo.h)**2 + data["M500"] *= (self._hdata / self.cosmo.h)**2 + # Get the 10s back in + data["L500"] *= 1e44 # ergs/s + data["M500"] *= 1e14 # Msun + + if max_redshift is not None: + data = data["z" <= max_redshift] + + self._data = data + + +class TwoMPPGalaxies(BaseSurvey): + """ + The 2M++ galaxy redshift catalogue [1], with the catalogue at [2]. + Removes fake galaxies used to fill the zone of avoidance. Note that the + stated redshift is in the CMB frame. + + Parameters + ---------- + fpath : str + File path to the catalogue. + + References + ---------- + [1] The 2M++ galaxy redshift catalogue; Lavaux, Guilhem, Hudson, Michael J. + [2] https://cdsarc.cds.unistra.fr/viz-bin/cat/J/MNRAS/416/2840#/article + [3] Improving NASA/IPAC Extragalactic Database Redshift Calculations + (2021); Anthony Carr and Tamara Davis + """ + + def __init__(self, fpath): + self._set_data(fpath) + + def _set_data(self, fpath): + """ + Set the catalogue + """ + from scipy.constants import c + # Read the catalogue and select non-fake galaxies + cat = numpy.genfromtxt(fpath, delimiter="|", ) + cat = cat[cat[:, 12] == 0, :] + # Pre=allocate array and fillt it + cols = [("RA", F64), ("DEC", F64), ("Ksmag", F64), ("ZCMB", F64), + ("DIST", F64)] + data = cols_to_structured(cat.shape[0], cols) + data["RA"] = cat[:, 1] + data["DEC"] = cat[:, 2] + data["Ksmag"] = cat[:, 5] + data["ZCMB"] = cat[:, 7] / (c * 1e-3) + self._data = data + + +class TwoMPPGroups(BaseSurvey): + """ + The 2M++ galaxy group catalogue [1], with the catalogue at [2]. + + Parameters + ---------- + fpath : str + File path to the catalogue. + + References + ---------- + [1] The 2M++ galaxy redshift catalogue; Lavaux, Guilhem, Hudson, Michael J. + [2] https://cdsarc.cds.unistra.fr/viz-bin/cat/J/MNRAS/416/2840#/article + [3] Improving NASA/IPAC Extragalactic Database Redshift Calculations + (2021); Anthony Carr and Tamara Davis + """ + + def __init__(self, fpath): + self._set_data(fpath) + + def _set_data(self, fpath): + """ + Set the catalogue + """ + cat = numpy.genfromtxt(fpath, delimiter="|", ) + # Pre-allocate and fill the array + cols = [("RA", F64), ("DEC", F64), ("K2mag", F64), + ("Rich", numpy.int64), ("sigma", F64)] + data = cols_to_structured(cat.shape[0], cols) + data["K2mag"] = cat[:, 3] + data["Rich"] = cat[:, 4] + data["sigma"] = cat[:, 7] + + # Convert galactic coordinates to RA, dec + glon = data[:, 1] + glat = data[:, 2] + coords = SkyCoord(l=glon*u.degree, b=glat*u.degree, frame='galactic') + coords = coords.transform_to("icrs") + data["RA"] = coords.ra + data["DEC"] = coords.dec + self._data = data diff --git a/csiborgtools/io/readsim.py b/csiborgtools/read/readsim.py similarity index 94% rename from csiborgtools/io/readsim.py rename to csiborgtools/read/readsim.py index 777b3ce..bb3a96f 100644 --- a/csiborgtools/io/readsim.py +++ b/csiborgtools/read/readsim.py @@ -22,8 +22,7 @@ from os import listdir from os.path import (join, isfile) from glob import glob from tqdm import tqdm - -from ..utils import (cols_to_structured, add_columns) +from ..utils import cols_to_structured F16 = numpy.float16 @@ -512,30 +511,3 @@ def read_mmain(n, srcdir, fname="Mmain_{}.npy"): out[name] = arr[:, i] return out - - -def merge_mmain_to_clumps(clumps, mmain): - """ - Merge columns from the `mmain` files to the `clump` file, matches them - by their halo index while assuming that the indices `index` in both arrays - are sorted. - - Parameters - ---------- - clumps : structured array - Clumps structured array. - mmain : structured array - Parent halo array whose information is to be merged into `clumps`. - - Returns - ------- - out : structured array - Array with added columns. - """ - X = numpy.full((clumps.size, 2), numpy.nan) - # Mask of which clumps have a mmain index - mask = numpy.isin(clumps["index"], mmain["index"]) - - X[mask, 0] = mmain["mass_cl"] - X[mask, 1] = mmain["sub_frac"] - return add_columns(clumps, X, ["mass_mmain", "sub_frac"]) diff --git a/csiborgtools/units/__init__.py b/csiborgtools/units/__init__.py index 0796f7b..9e4c738 100644 --- a/csiborgtools/units/__init__.py +++ b/csiborgtools/units/__init__.py @@ -14,4 +14,4 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. from .transforms import cartesian_to_radec # noqa -from .box_units import (BoxUnits, convert_from_boxunits) # noqa +from .box_units import (BoxUnits) # noqa diff --git a/csiborgtools/units/box_units.py b/csiborgtools/units/box_units.py index d090f22..9bd05d0 100644 --- a/csiborgtools/units/box_units.py +++ b/csiborgtools/units/box_units.py @@ -20,7 +20,7 @@ import numpy from scipy.interpolate import interp1d from astropy.cosmology import LambdaCDM from astropy import (constants, units) -from ..io import read_info +from ..read import read_info # Map of unit conversions @@ -329,65 +329,59 @@ class BoxUnits: return (density / self._unit_d * self._Msuncgs / (units.Mpc.to(units.cm))**3) + def convert_from_boxunits(self, data, names): + r""" + Convert columns named `names` in array `data` from box units to + physical units, such that + - length -> :math:`Mpc`, + - mass -> :math:`M_\odot`, + - density -> :math:`M_\odot / \mathrm{Mpc}^3`. -def convert_from_boxunits(data, names, boxunits): - r""" - Convert columns named `names` in array `data` from box units to physical - units, such that - - length -> :math:`Mpc`, - - mass -> :math:`M_\odot`, - - density -> :math:`M_\odot / \mathrm{Mpc}^3`. - Any other conversions are currently not implemented. Note that the array - is passed by reference and directly modified, even though it is also - explicitly returned. Additionally centres the box coordinates on the - observer, if they are being transformed. + Any other conversions are currently not implemented. Note that the + array is passed by reference and directly modified, even though it is + also explicitly returned. Additionally centres the box coordinates on + the observer, if they are being transformed. - Parameters - ---------- - data : structured array - Input array. - names : list of str - Columns to be converted. - boxunits : `BoxUnits` - Box units class of the simulation and snapshot. + Parameters + ---------- + data : structured array + Input array. + names : list of str + Columns to be converted. - Returns - ------- - data : structured array - Input array with converted columns. - """ - if not isinstance(boxunits, BoxUnits): - raise TypeError("`boxunits` must be of type `{}`. Currently `{}`." - .format(BoxUnits, type(boxunits))) - names = [names] if isinstance(names, str) else names + Returns + ------- + data : structured array + Input array with converted columns. + """ + names = [names] if isinstance(names, str) else names - # Shortcut for the transform functions - transforms = { - "length": boxunits.box2mpc, - "mass": boxunits.box2solarmass, - "density": boxunits.box2dens - } + # Shortcut for the transform functions + transforms = { + "length": self.box2mpc, + "mass": self.box2solarmass, + "density": self.box2dens + } - for name in names: - # Check that the name is even in the array - if name not in data.dtype.names: - raise ValueError("Name `{}` is not in `data` array.".format(name)) + for name in names: + # Check that the name is even in the array + if name not in data.dtype.names: + raise ValueError("Name `{}` not in `data` array.".format(name)) - # Convert - found = False - for unittype, suppnames in CONV_NAME.items(): - if name in suppnames: - data[name] = transforms[unittype](data[name]) - found = True - continue - # If nothing found - if not found: - raise NotImplementedError( - "Conversion of `{}` is not defined.".format(name)) + # Convert + found = False + for unittype, suppnames in CONV_NAME.items(): + if name in suppnames: + data[name] = transforms[unittype](data[name]) + found = True + continue + # If nothing found + if not found: + raise NotImplementedError( + "Conversion of `{}` is not defined.".format(name)) - # Center at the observer - if name in ["peak_x", "peak_y", "peak_z"]: - data[name] -= transforms["length"](0.5) - data[name] -= (0.5) + # Center at the observer + if name in ["peak_x", "peak_y", "peak_z"]: + data[name] -= transforms["length"](0.5) - return data + return data diff --git a/scripts/playground.ipynb b/scripts/playground.ipynb index bb35583..30abc0d 100644 --- a/scripts/playground.ipynb +++ b/scripts/playground.ipynb @@ -2,15 +2,24 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 160, "id": "5a38ed25", "metadata": { "ExecuteTime": { - "end_time": "2022-11-20T12:48:39.822802Z", - "start_time": "2022-11-20T12:48:03.992510Z" + "end_time": "2022-11-22T09:56:42.956157Z", + "start_time": "2022-11-22T09:56:41.999024Z" } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], "source": [ "import numpy as np\n", "%matplotlib notebook\n", @@ -26,12 +35,2603 @@ "%autoreload 2\n", "\n", "import joblib\n", - "from os.path import join" + "from os.path import join\n", + "\n", + "from scipy.interpolate import interp1d" ] }, + { + "cell_type": "code", + "execution_count": 163, + "id": "fa492543", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-22T09:57:21.730075Z", + "start_time": "2022-11-22T09:57:21.281306Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "'/mnt/extraspace/hdesmond/ramses_out_7636'" + ] + }, + "execution_count": 163, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sim = 7636\n", + "\n", + "csiborgtools.read.get_sim_path(sim)" + ] + }, + { + "cell_type": "code", + "execution_count": 162, + "id": "95759c3e", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-22T09:57:00.423419Z", + "start_time": "2022-11-22T09:57:00.337146Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "7636" + ] + }, + "execution_count": 162, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3c4c223", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "464b91b3", + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": 2, + "id": "2c98d50d", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T19:32:19.115273Z", + "start_time": "2022-11-21T19:29:12.079733Z" + }, + "scrolled": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [03:06<00:00, 18.69s/it]\n" + ] + } + ], + "source": [ + "cat = csiborgtools.read.CombinedHaloCatalogue()\n", + "# cat = csiborgtools.io.HaloCatalogue(7444, 951)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "72f406e5", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T19:35:07.774766Z", + "start_time": "2022-11-21T19:35:07.088097Z" + } + }, + "outputs": [], + "source": [ + "matcher = csiborgtools.match.RealisationsMatcher(cat)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "e165bef2", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T19:35:30.136090Z", + "start_time": "2022-11-21T19:35:09.690538Z" + }, + "scrolled": false + }, + "outputs": [], + "source": [ + "n = 0\n", + "match = matcher.cross_knn_position_single(n, 10)\n", + "# x.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 125, + "id": "47c34662", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T20:29:04.785499Z", + "start_time": "2022-11-21T20:29:04.046235Z" + } + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "/* global mpl */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function () {\n", + " if (typeof WebSocket !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof MozWebSocket !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert(\n", + " 'Your browser does not have WebSocket support. ' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.'\n", + " );\n", + " }\n", + "};\n", + "\n", + "mpl.figure = function (figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = this.ws.binaryType !== undefined;\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById('mpl-warnings');\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent =\n", + " 'This browser does not support binary websocket messages. ' +\n", + " 'Performance may be slow.';\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = document.createElement('div');\n", + " this.root.setAttribute('style', 'display: inline-block');\n", + " this._root_extra_style(this.root);\n", + "\n", + " parent_element.appendChild(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message('supports_binary', { value: fig.supports_binary });\n", + " fig.send_message('send_image_mode', {});\n", + " if (fig.ratio !== 1) {\n", + " fig.send_message('set_device_pixel_ratio', {\n", + " device_pixel_ratio: fig.ratio,\n", + " });\n", + " }\n", + " fig.send_message('refresh', {});\n", + " };\n", + "\n", + " this.imageObj.onload = function () {\n", + " if (fig.image_mode === 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function () {\n", + " fig.ws.close();\n", + " };\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "};\n", + "\n", + "mpl.figure.prototype._init_header = function () {\n", + " var titlebar = document.createElement('div');\n", + " titlebar.classList =\n", + " 'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n", + " var titletext = document.createElement('div');\n", + " titletext.classList = 'ui-dialog-title';\n", + " titletext.setAttribute(\n", + " 'style',\n", + " 'width: 100%; text-align: center; padding: 3px;'\n", + " );\n", + " titlebar.appendChild(titletext);\n", + " this.root.appendChild(titlebar);\n", + " this.header = titletext;\n", + "};\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n", + "\n", + "mpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n", + "\n", + "mpl.figure.prototype._init_canvas = function () {\n", + " var fig = this;\n", + "\n", + " var canvas_div = (this.canvas_div = document.createElement('div'));\n", + " canvas_div.setAttribute(\n", + " 'style',\n", + " 'border: 1px solid #ddd;' +\n", + " 'box-sizing: content-box;' +\n", + " 'clear: both;' +\n", + " 'min-height: 1px;' +\n", + " 'min-width: 1px;' +\n", + " 'outline: 0;' +\n", + " 'overflow: hidden;' +\n", + " 'position: relative;' +\n", + " 'resize: both;'\n", + " );\n", + "\n", + " function on_keyboard_event_closure(name) {\n", + " return function (event) {\n", + " return fig.key_event(event, name);\n", + " };\n", + " }\n", + "\n", + " canvas_div.addEventListener(\n", + " 'keydown',\n", + " on_keyboard_event_closure('key_press')\n", + " );\n", + " canvas_div.addEventListener(\n", + " 'keyup',\n", + " on_keyboard_event_closure('key_release')\n", + " );\n", + "\n", + " this._canvas_extra_style(canvas_div);\n", + " this.root.appendChild(canvas_div);\n", + "\n", + " var canvas = (this.canvas = document.createElement('canvas'));\n", + " canvas.classList.add('mpl-canvas');\n", + " canvas.setAttribute('style', 'box-sizing: content-box;');\n", + "\n", + " this.context = canvas.getContext('2d');\n", + "\n", + " var backingStore =\n", + " this.context.backingStorePixelRatio ||\n", + " this.context.webkitBackingStorePixelRatio ||\n", + " this.context.mozBackingStorePixelRatio ||\n", + " this.context.msBackingStorePixelRatio ||\n", + " this.context.oBackingStorePixelRatio ||\n", + " this.context.backingStorePixelRatio ||\n", + " 1;\n", + "\n", + " this.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband_canvas = (this.rubberband_canvas = document.createElement(\n", + " 'canvas'\n", + " ));\n", + " rubberband_canvas.setAttribute(\n", + " 'style',\n", + " 'box-sizing: content-box; position: absolute; left: 0; top: 0; z-index: 1;'\n", + " );\n", + "\n", + " // Apply a ponyfill if ResizeObserver is not implemented by browser.\n", + " if (this.ResizeObserver === undefined) {\n", + " if (window.ResizeObserver !== undefined) {\n", + " this.ResizeObserver = window.ResizeObserver;\n", + " } else {\n", + " var obs = _JSXTOOLS_RESIZE_OBSERVER({});\n", + " this.ResizeObserver = obs.ResizeObserver;\n", + " }\n", + " }\n", + "\n", + " this.resizeObserverInstance = new this.ResizeObserver(function (entries) {\n", + " var nentries = entries.length;\n", + " for (var i = 0; i < nentries; i++) {\n", + " var entry = entries[i];\n", + " var width, height;\n", + " if (entry.contentBoxSize) {\n", + " if (entry.contentBoxSize instanceof Array) {\n", + " // Chrome 84 implements new version of spec.\n", + " width = entry.contentBoxSize[0].inlineSize;\n", + " height = entry.contentBoxSize[0].blockSize;\n", + " } else {\n", + " // Firefox implements old version of spec.\n", + " width = entry.contentBoxSize.inlineSize;\n", + " height = entry.contentBoxSize.blockSize;\n", + " }\n", + " } else {\n", + " // Chrome <84 implements even older version of spec.\n", + " width = entry.contentRect.width;\n", + " height = entry.contentRect.height;\n", + " }\n", + "\n", + " // Keep the size of the canvas and rubber band canvas in sync with\n", + " // the canvas container.\n", + " if (entry.devicePixelContentBoxSize) {\n", + " // Chrome 84 implements new version of spec.\n", + " canvas.setAttribute(\n", + " 'width',\n", + " entry.devicePixelContentBoxSize[0].inlineSize\n", + " );\n", + " canvas.setAttribute(\n", + " 'height',\n", + " entry.devicePixelContentBoxSize[0].blockSize\n", + " );\n", + " } else {\n", + " canvas.setAttribute('width', width * fig.ratio);\n", + " canvas.setAttribute('height', height * fig.ratio);\n", + " }\n", + " canvas.setAttribute(\n", + " 'style',\n", + " 'width: ' + width + 'px; height: ' + height + 'px;'\n", + " );\n", + "\n", + " rubberband_canvas.setAttribute('width', width);\n", + " rubberband_canvas.setAttribute('height', height);\n", + "\n", + " // And update the size in Python. We ignore the initial 0/0 size\n", + " // that occurs as the element is placed into the DOM, which should\n", + " // otherwise not happen due to the minimum size styling.\n", + " if (fig.ws.readyState == 1 && width != 0 && height != 0) {\n", + " fig.request_resize(width, height);\n", + " }\n", + " }\n", + " });\n", + " this.resizeObserverInstance.observe(canvas_div);\n", + "\n", + " function on_mouse_event_closure(name) {\n", + " return function (event) {\n", + " return fig.mouse_event(event, name);\n", + " };\n", + " }\n", + "\n", + " rubberband_canvas.addEventListener(\n", + " 'mousedown',\n", + " on_mouse_event_closure('button_press')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseup',\n", + " on_mouse_event_closure('button_release')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'dblclick',\n", + " on_mouse_event_closure('dblclick')\n", + " );\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband_canvas.addEventListener(\n", + " 'mousemove',\n", + " on_mouse_event_closure('motion_notify')\n", + " );\n", + "\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseenter',\n", + " on_mouse_event_closure('figure_enter')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseleave',\n", + " on_mouse_event_closure('figure_leave')\n", + " );\n", + "\n", + " canvas_div.addEventListener('wheel', function (event) {\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " on_mouse_event_closure('scroll')(event);\n", + " });\n", + "\n", + " canvas_div.appendChild(canvas);\n", + " canvas_div.appendChild(rubberband_canvas);\n", + "\n", + " this.rubberband_context = rubberband_canvas.getContext('2d');\n", + " this.rubberband_context.strokeStyle = '#000000';\n", + "\n", + " this._resize_canvas = function (width, height, forward) {\n", + " if (forward) {\n", + " canvas_div.style.width = width + 'px';\n", + " canvas_div.style.height = height + 'px';\n", + " }\n", + " };\n", + "\n", + " // Disable right mouse context menu.\n", + " this.rubberband_canvas.addEventListener('contextmenu', function (_e) {\n", + " event.preventDefault();\n", + " return false;\n", + " });\n", + "\n", + " function set_focus() {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "};\n", + "\n", + "mpl.figure.prototype._init_toolbar = function () {\n", + " var fig = this;\n", + "\n", + " var toolbar = document.createElement('div');\n", + " toolbar.classList = 'mpl-toolbar';\n", + " this.root.appendChild(toolbar);\n", + "\n", + " function on_click_closure(name) {\n", + " return function (_event) {\n", + " return fig.toolbar_button_onclick(name);\n", + " };\n", + " }\n", + "\n", + " function on_mouseover_closure(tooltip) {\n", + " return function (event) {\n", + " if (!event.currentTarget.disabled) {\n", + " return fig.toolbar_button_onmouseover(tooltip);\n", + " }\n", + " };\n", + " }\n", + "\n", + " fig.buttons = {};\n", + " var buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'mpl-button-group';\n", + " for (var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " /* Instead of a spacer, we start a new button group. */\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + " buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'mpl-button-group';\n", + " continue;\n", + " }\n", + "\n", + " var button = (fig.buttons[name] = document.createElement('button'));\n", + " button.classList = 'mpl-widget';\n", + " button.setAttribute('role', 'button');\n", + " button.setAttribute('aria-disabled', 'false');\n", + " button.addEventListener('click', on_click_closure(method_name));\n", + " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", + "\n", + " var icon_img = document.createElement('img');\n", + " icon_img.src = '_images/' + image + '.png';\n", + " icon_img.srcset = '_images/' + image + '_large.png 2x';\n", + " icon_img.alt = tooltip;\n", + " button.appendChild(icon_img);\n", + "\n", + " buttonGroup.appendChild(button);\n", + " }\n", + "\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + "\n", + " var fmt_picker = document.createElement('select');\n", + " fmt_picker.classList = 'mpl-widget';\n", + " toolbar.appendChild(fmt_picker);\n", + " this.format_dropdown = fmt_picker;\n", + "\n", + " for (var ind in mpl.extensions) {\n", + " var fmt = mpl.extensions[ind];\n", + " var option = document.createElement('option');\n", + " option.selected = fmt === mpl.default_extension;\n", + " option.innerHTML = fmt;\n", + " fmt_picker.appendChild(option);\n", + " }\n", + "\n", + " var status_bar = document.createElement('span');\n", + " status_bar.classList = 'mpl-message';\n", + " toolbar.appendChild(status_bar);\n", + " this.message = status_bar;\n", + "};\n", + "\n", + "mpl.figure.prototype.request_resize = function (x_pixels, y_pixels) {\n", + " // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n", + " // which will in turn request a refresh of the image.\n", + " this.send_message('resize', { width: x_pixels, height: y_pixels });\n", + "};\n", + "\n", + "mpl.figure.prototype.send_message = function (type, properties) {\n", + " properties['type'] = type;\n", + " properties['figure_id'] = this.id;\n", + " this.ws.send(JSON.stringify(properties));\n", + "};\n", + "\n", + "mpl.figure.prototype.send_draw_message = function () {\n", + " if (!this.waiting) {\n", + " this.waiting = true;\n", + " this.ws.send(JSON.stringify({ type: 'draw', figure_id: this.id }));\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", + " var format_dropdown = fig.format_dropdown;\n", + " var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n", + " fig.ondownload(fig, format);\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_resize = function (fig, msg) {\n", + " var size = msg['size'];\n", + " if (size[0] !== fig.canvas.width || size[1] !== fig.canvas.height) {\n", + " fig._resize_canvas(size[0], size[1], msg['forward']);\n", + " fig.send_message('refresh', {});\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_rubberband = function (fig, msg) {\n", + " var x0 = msg['x0'] / fig.ratio;\n", + " var y0 = (fig.canvas.height - msg['y0']) / fig.ratio;\n", + " var x1 = msg['x1'] / fig.ratio;\n", + " var y1 = (fig.canvas.height - msg['y1']) / fig.ratio;\n", + " x0 = Math.floor(x0) + 0.5;\n", + " y0 = Math.floor(y0) + 0.5;\n", + " x1 = Math.floor(x1) + 0.5;\n", + " y1 = Math.floor(y1) + 0.5;\n", + " var min_x = Math.min(x0, x1);\n", + " var min_y = Math.min(y0, y1);\n", + " var width = Math.abs(x1 - x0);\n", + " var height = Math.abs(y1 - y0);\n", + "\n", + " fig.rubberband_context.clearRect(\n", + " 0,\n", + " 0,\n", + " fig.canvas.width / fig.ratio,\n", + " fig.canvas.height / fig.ratio\n", + " );\n", + "\n", + " fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_figure_label = function (fig, msg) {\n", + " // Updates the figure title.\n", + " fig.header.textContent = msg['label'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_cursor = function (fig, msg) {\n", + " fig.rubberband_canvas.style.cursor = msg['cursor'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_message = function (fig, msg) {\n", + " fig.message.textContent = msg['message'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_draw = function (fig, _msg) {\n", + " // Request the server to send over a new figure.\n", + " fig.send_draw_message();\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_image_mode = function (fig, msg) {\n", + " fig.image_mode = msg['mode'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_history_buttons = function (fig, msg) {\n", + " for (var key in msg) {\n", + " if (!(key in fig.buttons)) {\n", + " continue;\n", + " }\n", + " fig.buttons[key].disabled = !msg[key];\n", + " fig.buttons[key].setAttribute('aria-disabled', !msg[key]);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_navigate_mode = function (fig, msg) {\n", + " if (msg['mode'] === 'PAN') {\n", + " fig.buttons['Pan'].classList.add('active');\n", + " fig.buttons['Zoom'].classList.remove('active');\n", + " } else if (msg['mode'] === 'ZOOM') {\n", + " fig.buttons['Pan'].classList.remove('active');\n", + " fig.buttons['Zoom'].classList.add('active');\n", + " } else {\n", + " fig.buttons['Pan'].classList.remove('active');\n", + " fig.buttons['Zoom'].classList.remove('active');\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.updated_canvas_event = function () {\n", + " // Called whenever the canvas gets updated.\n", + " this.send_message('ack', {});\n", + "};\n", + "\n", + "// A function to construct a web socket function for onmessage handling.\n", + "// Called in the figure constructor.\n", + "mpl.figure.prototype._make_on_message_function = function (fig) {\n", + " return function socket_on_message(evt) {\n", + " if (evt.data instanceof Blob) {\n", + " var img = evt.data;\n", + " if (img.type !== 'image/png') {\n", + " /* FIXME: We get \"Resource interpreted as Image but\n", + " * transferred with MIME type text/plain:\" errors on\n", + " * Chrome. But how to set the MIME type? It doesn't seem\n", + " * to be part of the websocket stream */\n", + " img.type = 'image/png';\n", + " }\n", + "\n", + " /* Free the memory for the previous frames */\n", + " if (fig.imageObj.src) {\n", + " (window.URL || window.webkitURL).revokeObjectURL(\n", + " fig.imageObj.src\n", + " );\n", + " }\n", + "\n", + " fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n", + " img\n", + " );\n", + " fig.updated_canvas_event();\n", + " fig.waiting = false;\n", + " return;\n", + " } else if (\n", + " typeof evt.data === 'string' &&\n", + " evt.data.slice(0, 21) === 'data:image/png;base64'\n", + " ) {\n", + " fig.imageObj.src = evt.data;\n", + " fig.updated_canvas_event();\n", + " fig.waiting = false;\n", + " return;\n", + " }\n", + "\n", + " var msg = JSON.parse(evt.data);\n", + " var msg_type = msg['type'];\n", + "\n", + " // Call the \"handle_{type}\" callback, which takes\n", + " // the figure and JSON message as its only arguments.\n", + " try {\n", + " var callback = fig['handle_' + msg_type];\n", + " } catch (e) {\n", + " console.log(\n", + " \"No handler for the '\" + msg_type + \"' message type: \",\n", + " msg\n", + " );\n", + " return;\n", + " }\n", + "\n", + " if (callback) {\n", + " try {\n", + " // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n", + " callback(fig, msg);\n", + " } catch (e) {\n", + " console.log(\n", + " \"Exception inside the 'handler_\" + msg_type + \"' callback:\",\n", + " e,\n", + " e.stack,\n", + " msg\n", + " );\n", + " }\n", + " }\n", + " };\n", + "};\n", + "\n", + "// from https://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n", + "mpl.findpos = function (e) {\n", + " //this section is from http://www.quirksmode.org/js/events_properties.html\n", + " var targ;\n", + " if (!e) {\n", + " e = window.event;\n", + " }\n", + " if (e.target) {\n", + " targ = e.target;\n", + " } else if (e.srcElement) {\n", + " targ = e.srcElement;\n", + " }\n", + " if (targ.nodeType === 3) {\n", + " // defeat Safari bug\n", + " targ = targ.parentNode;\n", + " }\n", + "\n", + " // pageX,Y are the mouse positions relative to the document\n", + " var boundingRect = targ.getBoundingClientRect();\n", + " var x = e.pageX - (boundingRect.left + document.body.scrollLeft);\n", + " var y = e.pageY - (boundingRect.top + document.body.scrollTop);\n", + "\n", + " return { x: x, y: y };\n", + "};\n", + "\n", + "/*\n", + " * return a copy of an object with only non-object keys\n", + " * we need this to avoid circular references\n", + " * https://stackoverflow.com/a/24161582/3208463\n", + " */\n", + "function simpleKeys(original) {\n", + " return Object.keys(original).reduce(function (obj, key) {\n", + " if (typeof original[key] !== 'object') {\n", + " obj[key] = original[key];\n", + " }\n", + " return obj;\n", + " }, {});\n", + "}\n", + "\n", + "mpl.figure.prototype.mouse_event = function (event, name) {\n", + " var canvas_pos = mpl.findpos(event);\n", + "\n", + " if (name === 'button_press') {\n", + " this.canvas.focus();\n", + " this.canvas_div.focus();\n", + " }\n", + "\n", + " var x = canvas_pos.x * this.ratio;\n", + " var y = canvas_pos.y * this.ratio;\n", + "\n", + " this.send_message(name, {\n", + " x: x,\n", + " y: y,\n", + " button: event.button,\n", + " step: event.step,\n", + " guiEvent: simpleKeys(event),\n", + " });\n", + "\n", + " /* This prevents the web browser from automatically changing to\n", + " * the text insertion cursor when the button is pressed. We want\n", + " * to control all of the cursor setting manually through the\n", + " * 'cursor' event from matplotlib */\n", + " event.preventDefault();\n", + " return false;\n", + "};\n", + "\n", + "mpl.figure.prototype._key_event_extra = function (_event, _name) {\n", + " // Handle any extra behaviour associated with a key event\n", + "};\n", + "\n", + "mpl.figure.prototype.key_event = function (event, name) {\n", + " // Prevent repeat events\n", + " if (name === 'key_press') {\n", + " if (event.key === this._key) {\n", + " return;\n", + " } else {\n", + " this._key = event.key;\n", + " }\n", + " }\n", + " if (name === 'key_release') {\n", + " this._key = null;\n", + " }\n", + "\n", + " var value = '';\n", + " if (event.ctrlKey && event.key !== 'Control') {\n", + " value += 'ctrl+';\n", + " }\n", + " else if (event.altKey && event.key !== 'Alt') {\n", + " value += 'alt+';\n", + " }\n", + " else if (event.shiftKey && event.key !== 'Shift') {\n", + " value += 'shift+';\n", + " }\n", + "\n", + " value += 'k' + event.key;\n", + "\n", + " this._key_event_extra(event, name);\n", + "\n", + " this.send_message(name, { key: value, guiEvent: simpleKeys(event) });\n", + " return false;\n", + "};\n", + "\n", + "mpl.figure.prototype.toolbar_button_onclick = function (name) {\n", + " if (name === 'download') {\n", + " this.handle_save(this, null);\n", + " } else {\n", + " this.send_message('toolbar_button', { name: name });\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.toolbar_button_onmouseover = function (tooltip) {\n", + " this.message.textContent = tooltip;\n", + "};\n", + "\n", + "///////////////// REMAINING CONTENT GENERATED BY embed_js.py /////////////////\n", + "// prettier-ignore\n", + "var _JSXTOOLS_RESIZE_OBSERVER=function(A){var t,i=new WeakMap,n=new WeakMap,a=new WeakMap,r=new WeakMap,o=new Set;function s(e){if(!(this instanceof s))throw new TypeError(\"Constructor requires 'new' operator\");i.set(this,e)}function h(){throw new TypeError(\"Function is not a constructor\")}function c(e,t,i,n){e=0 in arguments?Number(arguments[0]):0,t=1 in arguments?Number(arguments[1]):0,i=2 in arguments?Number(arguments[2]):0,n=3 in arguments?Number(arguments[3]):0,this.right=(this.x=this.left=e)+(this.width=i),this.bottom=(this.y=this.top=t)+(this.height=n),Object.freeze(this)}function d(){t=requestAnimationFrame(d);var s=new WeakMap,p=new Set;o.forEach((function(t){r.get(t).forEach((function(i){var r=t instanceof window.SVGElement,o=a.get(t),d=r?0:parseFloat(o.paddingTop),f=r?0:parseFloat(o.paddingRight),l=r?0:parseFloat(o.paddingBottom),u=r?0:parseFloat(o.paddingLeft),g=r?0:parseFloat(o.borderTopWidth),m=r?0:parseFloat(o.borderRightWidth),w=r?0:parseFloat(o.borderBottomWidth),b=u+f,F=d+l,v=(r?0:parseFloat(o.borderLeftWidth))+m,W=g+w,y=r?0:t.offsetHeight-W-t.clientHeight,E=r?0:t.offsetWidth-v-t.clientWidth,R=b+v,z=F+W,M=r?t.width:parseFloat(o.width)-R-E,O=r?t.height:parseFloat(o.height)-z-y;if(n.has(t)){var k=n.get(t);if(k[0]===M&&k[1]===O)return}n.set(t,[M,O]);var S=Object.create(h.prototype);S.target=t,S.contentRect=new c(u,d,M,O),s.has(i)||(s.set(i,[]),p.add(i)),s.get(i).push(S)}))})),p.forEach((function(e){i.get(e).call(e,s.get(e),e)}))}return s.prototype.observe=function(i){if(i instanceof window.Element){r.has(i)||(r.set(i,new Set),o.add(i),a.set(i,window.getComputedStyle(i)));var n=r.get(i);n.has(this)||n.add(this),cancelAnimationFrame(t),t=requestAnimationFrame(d)}},s.prototype.unobserve=function(i){if(i instanceof window.Element&&r.has(i)){var n=r.get(i);n.has(this)&&(n.delete(this),n.size||(r.delete(i),o.delete(i))),n.size||r.delete(i),o.size||cancelAnimationFrame(t)}},A.DOMRectReadOnly=c,A.ResizeObserver=s,A.ResizeObserverEntry=h,A}; // eslint-disable-line\n", + "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Left button pans, Right button zooms\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-arrows\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\\nx/y fixes axis\", \"fa fa-square-o\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o\", \"download\"]];\n", + "\n", + "mpl.extensions = [\"eps\", \"jpeg\", \"pgf\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\", \"webp\"];\n", + "\n", + "mpl.default_extension = \"png\";/* global mpl */\n", + "\n", + "var comm_websocket_adapter = function (comm) {\n", + " // Create a \"websocket\"-like object which calls the given IPython comm\n", + " // object with the appropriate methods. Currently this is a non binary\n", + " // socket, so there is still some room for performance tuning.\n", + " var ws = {};\n", + "\n", + " ws.binaryType = comm.kernel.ws.binaryType;\n", + " ws.readyState = comm.kernel.ws.readyState;\n", + " function updateReadyState(_event) {\n", + " if (comm.kernel.ws) {\n", + " ws.readyState = comm.kernel.ws.readyState;\n", + " } else {\n", + " ws.readyState = 3; // Closed state.\n", + " }\n", + " }\n", + " comm.kernel.ws.addEventListener('open', updateReadyState);\n", + " comm.kernel.ws.addEventListener('close', updateReadyState);\n", + " comm.kernel.ws.addEventListener('error', updateReadyState);\n", + "\n", + " ws.close = function () {\n", + " comm.close();\n", + " };\n", + " ws.send = function (m) {\n", + " //console.log('sending', m);\n", + " comm.send(m);\n", + " };\n", + " // Register the callback with on_msg.\n", + " comm.on_msg(function (msg) {\n", + " //console.log('receiving', msg['content']['data'], msg);\n", + " var data = msg['content']['data'];\n", + " if (data['blob'] !== undefined) {\n", + " data = {\n", + " data: new Blob(msg['buffers'], { type: data['blob'] }),\n", + " };\n", + " }\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", + " ws.onmessage(data);\n", + " });\n", + " return ws;\n", + "};\n", + "\n", + "mpl.mpl_figure_comm = function (comm, msg) {\n", + " // This is the function which gets called when the mpl process\n", + " // starts-up an IPython Comm through the \"matplotlib\" channel.\n", + "\n", + " var id = msg.content.data.id;\n", + " // Get hold of the div created by the display call when the Comm\n", + " // socket was opened in Python.\n", + " var element = document.getElementById(id);\n", + " var ws_proxy = comm_websocket_adapter(comm);\n", + "\n", + " function ondownload(figure, _format) {\n", + " window.open(figure.canvas.toDataURL());\n", + " }\n", + "\n", + " var fig = new mpl.figure(id, ws_proxy, ondownload, element);\n", + "\n", + " // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n", + " // web socket which is closed, not our websocket->open comm proxy.\n", + " ws_proxy.onopen();\n", + "\n", + " fig.parent_element = element;\n", + " fig.cell_info = mpl.find_output_cell(\"
\");\n", + " if (!fig.cell_info) {\n", + " console.error('Failed to find cell for figure', id, fig);\n", + " return;\n", + " }\n", + " fig.cell_info[0].output_area.element.on(\n", + " 'cleared',\n", + " { fig: fig },\n", + " fig._remove_fig_handler\n", + " );\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_close = function (fig, msg) {\n", + " var width = fig.canvas.width / fig.ratio;\n", + " fig.cell_info[0].output_area.element.off(\n", + " 'cleared',\n", + " fig._remove_fig_handler\n", + " );\n", + " fig.resizeObserverInstance.unobserve(fig.canvas_div);\n", + "\n", + " // Update the output cell to use the data from the current canvas.\n", + " fig.push_to_output();\n", + " var dataURL = fig.canvas.toDataURL();\n", + " // Re-enable the keyboard manager in IPython - without this line, in FF,\n", + " // the notebook keyboard shortcuts fail.\n", + " IPython.keyboard_manager.enable();\n", + " fig.parent_element.innerHTML =\n", + " '';\n", + " fig.close_ws(fig, msg);\n", + "};\n", + "\n", + "mpl.figure.prototype.close_ws = function (fig, msg) {\n", + " fig.send_message('closing', msg);\n", + " // fig.ws.close()\n", + "};\n", + "\n", + "mpl.figure.prototype.push_to_output = function (_remove_interactive) {\n", + " // Turn the data on the canvas into data in the output cell.\n", + " var width = this.canvas.width / this.ratio;\n", + " var dataURL = this.canvas.toDataURL();\n", + " this.cell_info[1]['text/html'] =\n", + " '';\n", + "};\n", + "\n", + "mpl.figure.prototype.updated_canvas_event = function () {\n", + " // Tell IPython that the notebook contents must change.\n", + " IPython.notebook.set_dirty(true);\n", + " this.send_message('ack', {});\n", + " var fig = this;\n", + " // Wait a second, then push the new image to the DOM so\n", + " // that it is saved nicely (might be nice to debounce this).\n", + " setTimeout(function () {\n", + " fig.push_to_output();\n", + " }, 1000);\n", + "};\n", + "\n", + "mpl.figure.prototype._init_toolbar = function () {\n", + " var fig = this;\n", + "\n", + " var toolbar = document.createElement('div');\n", + " toolbar.classList = 'btn-toolbar';\n", + " this.root.appendChild(toolbar);\n", + "\n", + " function on_click_closure(name) {\n", + " return function (_event) {\n", + " return fig.toolbar_button_onclick(name);\n", + " };\n", + " }\n", + "\n", + " function on_mouseover_closure(tooltip) {\n", + " return function (event) {\n", + " if (!event.currentTarget.disabled) {\n", + " return fig.toolbar_button_onmouseover(tooltip);\n", + " }\n", + " };\n", + " }\n", + "\n", + " fig.buttons = {};\n", + " var buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'btn-group';\n", + " var button;\n", + " for (var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " /* Instead of a spacer, we start a new button group. */\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + " buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'btn-group';\n", + " continue;\n", + " }\n", + "\n", + " button = fig.buttons[name] = document.createElement('button');\n", + " button.classList = 'btn btn-default';\n", + " button.href = '#';\n", + " button.title = name;\n", + " button.innerHTML = '';\n", + " button.addEventListener('click', on_click_closure(method_name));\n", + " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", + " buttonGroup.appendChild(button);\n", + " }\n", + "\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = document.createElement('span');\n", + " status_bar.classList = 'mpl-message pull-right';\n", + " toolbar.appendChild(status_bar);\n", + " this.message = status_bar;\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = document.createElement('div');\n", + " buttongrp.classList = 'btn-group inline pull-right';\n", + " button = document.createElement('button');\n", + " button.classList = 'btn btn-mini btn-primary';\n", + " button.href = '#';\n", + " button.title = 'Stop Interaction';\n", + " button.innerHTML = '';\n", + " button.addEventListener('click', function (_evt) {\n", + " fig.handle_close(fig, {});\n", + " });\n", + " button.addEventListener(\n", + " 'mouseover',\n", + " on_mouseover_closure('Stop Interaction')\n", + " );\n", + " buttongrp.appendChild(button);\n", + " var titlebar = this.root.querySelector('.ui-dialog-titlebar');\n", + " titlebar.insertBefore(buttongrp, titlebar.firstChild);\n", + "};\n", + "\n", + "mpl.figure.prototype._remove_fig_handler = function (event) {\n", + " var fig = event.data.fig;\n", + " if (event.target !== this) {\n", + " // Ignore bubbled events from children.\n", + " return;\n", + " }\n", + " fig.close_ws(fig, {});\n", + "};\n", + "\n", + "mpl.figure.prototype._root_extra_style = function (el) {\n", + " el.style.boxSizing = 'content-box'; // override notebook setting of border-box.\n", + "};\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function (el) {\n", + " // this is important to make the div 'focusable\n", + " el.setAttribute('tabindex', 0);\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " } else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype._key_event_extra = function (event, _name) {\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which === 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", + " fig.ondownload(fig, null);\n", + "};\n", + "\n", + "mpl.find_output_cell = function (html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i = 0; i < ncells; i++) {\n", + " var cell = cells[i];\n", + " if (cell.cell_type === 'code') {\n", + " for (var j = 0; j < cell.output_area.outputs.length; j++) {\n", + " var data = cell.output_area.outputs[j];\n", + " if (data.data) {\n", + " // IPython >= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] === html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "};\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel !== null) {\n", + " IPython.notebook.kernel.comm_manager.register_target(\n", + " 'matplotlib',\n", + " mpl.mpl_figure_comm\n", + " );\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "q = np.linspace(0, 100, 10000)\n", + "ps = np.percentile(cat[n][\"m200\"], q)\n", + "f = interp1d(q, np.log10(ps))\n", + "\n", + "pcuts = [99.99, 95., 65.]\n", + "\n", + "\n", + "colors = plt.rcParams['axes.prop_cycle'].by_key()['color']\n", + "plt.figure()\n", + "plt.plot(ps, q)\n", + "\n", + "for i, pcut in enumerate(pcuts):\n", + " N = int((100 - pcut) / 100 * cat[n][\"m200\"].size)\n", + " plt.axvline(10**float(f(pcut)), ls=\"--\", c=colors[i + 1],\n", + " label=r\"${}^{{\\rm th}} percentile, {}$\".format(pcut, N), lw=1)\n", + " plt.axhline(pcut, c=colors[i + 1], ls=\"--\", lw=1)\n", + "plt.xscale(\"log\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 157, + "id": "de2a2714", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T20:30:43.637731Z", + "start_time": "2022-11-21T20:30:43.589062Z" + } + }, + "outputs": [], + "source": [ + "ids = np.where(cat[n][\"m200\"] > 10**f(pcuts[2]))[0]\n", + "ids = np.random.choice(ids, size=3)" + ] + }, + { + "cell_type": "code", + "execution_count": 159, + "id": "2ab2af25", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T20:30:54.506179Z", + "start_time": "2022-11-21T20:30:53.801870Z" + } + }, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "/* global mpl */\n", + "window.mpl = {};\n", + "\n", + "mpl.get_websocket_type = function () {\n", + " if (typeof WebSocket !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof MozWebSocket !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert(\n", + " 'Your browser does not have WebSocket support. ' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.'\n", + " );\n", + " }\n", + "};\n", + "\n", + "mpl.figure = function (figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = this.ws.binaryType !== undefined;\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById('mpl-warnings');\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent =\n", + " 'This browser does not support binary websocket messages. ' +\n", + " 'Performance may be slow.';\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = document.createElement('div');\n", + " this.root.setAttribute('style', 'display: inline-block');\n", + " this._root_extra_style(this.root);\n", + "\n", + " parent_element.appendChild(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message('supports_binary', { value: fig.supports_binary });\n", + " fig.send_message('send_image_mode', {});\n", + " if (fig.ratio !== 1) {\n", + " fig.send_message('set_device_pixel_ratio', {\n", + " device_pixel_ratio: fig.ratio,\n", + " });\n", + " }\n", + " fig.send_message('refresh', {});\n", + " };\n", + "\n", + " this.imageObj.onload = function () {\n", + " if (fig.image_mode === 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function () {\n", + " fig.ws.close();\n", + " };\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "};\n", + "\n", + "mpl.figure.prototype._init_header = function () {\n", + " var titlebar = document.createElement('div');\n", + " titlebar.classList =\n", + " 'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n", + " var titletext = document.createElement('div');\n", + " titletext.classList = 'ui-dialog-title';\n", + " titletext.setAttribute(\n", + " 'style',\n", + " 'width: 100%; text-align: center; padding: 3px;'\n", + " );\n", + " titlebar.appendChild(titletext);\n", + " this.root.appendChild(titlebar);\n", + " this.header = titletext;\n", + "};\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n", + "\n", + "mpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n", + "\n", + "mpl.figure.prototype._init_canvas = function () {\n", + " var fig = this;\n", + "\n", + " var canvas_div = (this.canvas_div = document.createElement('div'));\n", + " canvas_div.setAttribute(\n", + " 'style',\n", + " 'border: 1px solid #ddd;' +\n", + " 'box-sizing: content-box;' +\n", + " 'clear: both;' +\n", + " 'min-height: 1px;' +\n", + " 'min-width: 1px;' +\n", + " 'outline: 0;' +\n", + " 'overflow: hidden;' +\n", + " 'position: relative;' +\n", + " 'resize: both;'\n", + " );\n", + "\n", + " function on_keyboard_event_closure(name) {\n", + " return function (event) {\n", + " return fig.key_event(event, name);\n", + " };\n", + " }\n", + "\n", + " canvas_div.addEventListener(\n", + " 'keydown',\n", + " on_keyboard_event_closure('key_press')\n", + " );\n", + " canvas_div.addEventListener(\n", + " 'keyup',\n", + " on_keyboard_event_closure('key_release')\n", + " );\n", + "\n", + " this._canvas_extra_style(canvas_div);\n", + " this.root.appendChild(canvas_div);\n", + "\n", + " var canvas = (this.canvas = document.createElement('canvas'));\n", + " canvas.classList.add('mpl-canvas');\n", + " canvas.setAttribute('style', 'box-sizing: content-box;');\n", + "\n", + " this.context = canvas.getContext('2d');\n", + "\n", + " var backingStore =\n", + " this.context.backingStorePixelRatio ||\n", + " this.context.webkitBackingStorePixelRatio ||\n", + " this.context.mozBackingStorePixelRatio ||\n", + " this.context.msBackingStorePixelRatio ||\n", + " this.context.oBackingStorePixelRatio ||\n", + " this.context.backingStorePixelRatio ||\n", + " 1;\n", + "\n", + " this.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband_canvas = (this.rubberband_canvas = document.createElement(\n", + " 'canvas'\n", + " ));\n", + " rubberband_canvas.setAttribute(\n", + " 'style',\n", + " 'box-sizing: content-box; position: absolute; left: 0; top: 0; z-index: 1;'\n", + " );\n", + "\n", + " // Apply a ponyfill if ResizeObserver is not implemented by browser.\n", + " if (this.ResizeObserver === undefined) {\n", + " if (window.ResizeObserver !== undefined) {\n", + " this.ResizeObserver = window.ResizeObserver;\n", + " } else {\n", + " var obs = _JSXTOOLS_RESIZE_OBSERVER({});\n", + " this.ResizeObserver = obs.ResizeObserver;\n", + " }\n", + " }\n", + "\n", + " this.resizeObserverInstance = new this.ResizeObserver(function (entries) {\n", + " var nentries = entries.length;\n", + " for (var i = 0; i < nentries; i++) {\n", + " var entry = entries[i];\n", + " var width, height;\n", + " if (entry.contentBoxSize) {\n", + " if (entry.contentBoxSize instanceof Array) {\n", + " // Chrome 84 implements new version of spec.\n", + " width = entry.contentBoxSize[0].inlineSize;\n", + " height = entry.contentBoxSize[0].blockSize;\n", + " } else {\n", + " // Firefox implements old version of spec.\n", + " width = entry.contentBoxSize.inlineSize;\n", + " height = entry.contentBoxSize.blockSize;\n", + " }\n", + " } else {\n", + " // Chrome <84 implements even older version of spec.\n", + " width = entry.contentRect.width;\n", + " height = entry.contentRect.height;\n", + " }\n", + "\n", + " // Keep the size of the canvas and rubber band canvas in sync with\n", + " // the canvas container.\n", + " if (entry.devicePixelContentBoxSize) {\n", + " // Chrome 84 implements new version of spec.\n", + " canvas.setAttribute(\n", + " 'width',\n", + " entry.devicePixelContentBoxSize[0].inlineSize\n", + " );\n", + " canvas.setAttribute(\n", + " 'height',\n", + " entry.devicePixelContentBoxSize[0].blockSize\n", + " );\n", + " } else {\n", + " canvas.setAttribute('width', width * fig.ratio);\n", + " canvas.setAttribute('height', height * fig.ratio);\n", + " }\n", + " canvas.setAttribute(\n", + " 'style',\n", + " 'width: ' + width + 'px; height: ' + height + 'px;'\n", + " );\n", + "\n", + " rubberband_canvas.setAttribute('width', width);\n", + " rubberband_canvas.setAttribute('height', height);\n", + "\n", + " // And update the size in Python. We ignore the initial 0/0 size\n", + " // that occurs as the element is placed into the DOM, which should\n", + " // otherwise not happen due to the minimum size styling.\n", + " if (fig.ws.readyState == 1 && width != 0 && height != 0) {\n", + " fig.request_resize(width, height);\n", + " }\n", + " }\n", + " });\n", + " this.resizeObserverInstance.observe(canvas_div);\n", + "\n", + " function on_mouse_event_closure(name) {\n", + " return function (event) {\n", + " return fig.mouse_event(event, name);\n", + " };\n", + " }\n", + "\n", + " rubberband_canvas.addEventListener(\n", + " 'mousedown',\n", + " on_mouse_event_closure('button_press')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseup',\n", + " on_mouse_event_closure('button_release')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'dblclick',\n", + " on_mouse_event_closure('dblclick')\n", + " );\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband_canvas.addEventListener(\n", + " 'mousemove',\n", + " on_mouse_event_closure('motion_notify')\n", + " );\n", + "\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseenter',\n", + " on_mouse_event_closure('figure_enter')\n", + " );\n", + " rubberband_canvas.addEventListener(\n", + " 'mouseleave',\n", + " on_mouse_event_closure('figure_leave')\n", + " );\n", + "\n", + " canvas_div.addEventListener('wheel', function (event) {\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " on_mouse_event_closure('scroll')(event);\n", + " });\n", + "\n", + " canvas_div.appendChild(canvas);\n", + " canvas_div.appendChild(rubberband_canvas);\n", + "\n", + " this.rubberband_context = rubberband_canvas.getContext('2d');\n", + " this.rubberband_context.strokeStyle = '#000000';\n", + "\n", + " this._resize_canvas = function (width, height, forward) {\n", + " if (forward) {\n", + " canvas_div.style.width = width + 'px';\n", + " canvas_div.style.height = height + 'px';\n", + " }\n", + " };\n", + "\n", + " // Disable right mouse context menu.\n", + " this.rubberband_canvas.addEventListener('contextmenu', function (_e) {\n", + " event.preventDefault();\n", + " return false;\n", + " });\n", + "\n", + " function set_focus() {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "};\n", + "\n", + "mpl.figure.prototype._init_toolbar = function () {\n", + " var fig = this;\n", + "\n", + " var toolbar = document.createElement('div');\n", + " toolbar.classList = 'mpl-toolbar';\n", + " this.root.appendChild(toolbar);\n", + "\n", + " function on_click_closure(name) {\n", + " return function (_event) {\n", + " return fig.toolbar_button_onclick(name);\n", + " };\n", + " }\n", + "\n", + " function on_mouseover_closure(tooltip) {\n", + " return function (event) {\n", + " if (!event.currentTarget.disabled) {\n", + " return fig.toolbar_button_onmouseover(tooltip);\n", + " }\n", + " };\n", + " }\n", + "\n", + " fig.buttons = {};\n", + " var buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'mpl-button-group';\n", + " for (var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " /* Instead of a spacer, we start a new button group. */\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + " buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'mpl-button-group';\n", + " continue;\n", + " }\n", + "\n", + " var button = (fig.buttons[name] = document.createElement('button'));\n", + " button.classList = 'mpl-widget';\n", + " button.setAttribute('role', 'button');\n", + " button.setAttribute('aria-disabled', 'false');\n", + " button.addEventListener('click', on_click_closure(method_name));\n", + " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", + "\n", + " var icon_img = document.createElement('img');\n", + " icon_img.src = '_images/' + image + '.png';\n", + " icon_img.srcset = '_images/' + image + '_large.png 2x';\n", + " icon_img.alt = tooltip;\n", + " button.appendChild(icon_img);\n", + "\n", + " buttonGroup.appendChild(button);\n", + " }\n", + "\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + "\n", + " var fmt_picker = document.createElement('select');\n", + " fmt_picker.classList = 'mpl-widget';\n", + " toolbar.appendChild(fmt_picker);\n", + " this.format_dropdown = fmt_picker;\n", + "\n", + " for (var ind in mpl.extensions) {\n", + " var fmt = mpl.extensions[ind];\n", + " var option = document.createElement('option');\n", + " option.selected = fmt === mpl.default_extension;\n", + " option.innerHTML = fmt;\n", + " fmt_picker.appendChild(option);\n", + " }\n", + "\n", + " var status_bar = document.createElement('span');\n", + " status_bar.classList = 'mpl-message';\n", + " toolbar.appendChild(status_bar);\n", + " this.message = status_bar;\n", + "};\n", + "\n", + "mpl.figure.prototype.request_resize = function (x_pixels, y_pixels) {\n", + " // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n", + " // which will in turn request a refresh of the image.\n", + " this.send_message('resize', { width: x_pixels, height: y_pixels });\n", + "};\n", + "\n", + "mpl.figure.prototype.send_message = function (type, properties) {\n", + " properties['type'] = type;\n", + " properties['figure_id'] = this.id;\n", + " this.ws.send(JSON.stringify(properties));\n", + "};\n", + "\n", + "mpl.figure.prototype.send_draw_message = function () {\n", + " if (!this.waiting) {\n", + " this.waiting = true;\n", + " this.ws.send(JSON.stringify({ type: 'draw', figure_id: this.id }));\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", + " var format_dropdown = fig.format_dropdown;\n", + " var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n", + " fig.ondownload(fig, format);\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_resize = function (fig, msg) {\n", + " var size = msg['size'];\n", + " if (size[0] !== fig.canvas.width || size[1] !== fig.canvas.height) {\n", + " fig._resize_canvas(size[0], size[1], msg['forward']);\n", + " fig.send_message('refresh', {});\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_rubberband = function (fig, msg) {\n", + " var x0 = msg['x0'] / fig.ratio;\n", + " var y0 = (fig.canvas.height - msg['y0']) / fig.ratio;\n", + " var x1 = msg['x1'] / fig.ratio;\n", + " var y1 = (fig.canvas.height - msg['y1']) / fig.ratio;\n", + " x0 = Math.floor(x0) + 0.5;\n", + " y0 = Math.floor(y0) + 0.5;\n", + " x1 = Math.floor(x1) + 0.5;\n", + " y1 = Math.floor(y1) + 0.5;\n", + " var min_x = Math.min(x0, x1);\n", + " var min_y = Math.min(y0, y1);\n", + " var width = Math.abs(x1 - x0);\n", + " var height = Math.abs(y1 - y0);\n", + "\n", + " fig.rubberband_context.clearRect(\n", + " 0,\n", + " 0,\n", + " fig.canvas.width / fig.ratio,\n", + " fig.canvas.height / fig.ratio\n", + " );\n", + "\n", + " fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_figure_label = function (fig, msg) {\n", + " // Updates the figure title.\n", + " fig.header.textContent = msg['label'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_cursor = function (fig, msg) {\n", + " fig.rubberband_canvas.style.cursor = msg['cursor'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_message = function (fig, msg) {\n", + " fig.message.textContent = msg['message'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_draw = function (fig, _msg) {\n", + " // Request the server to send over a new figure.\n", + " fig.send_draw_message();\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_image_mode = function (fig, msg) {\n", + " fig.image_mode = msg['mode'];\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_history_buttons = function (fig, msg) {\n", + " for (var key in msg) {\n", + " if (!(key in fig.buttons)) {\n", + " continue;\n", + " }\n", + " fig.buttons[key].disabled = !msg[key];\n", + " fig.buttons[key].setAttribute('aria-disabled', !msg[key]);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_navigate_mode = function (fig, msg) {\n", + " if (msg['mode'] === 'PAN') {\n", + " fig.buttons['Pan'].classList.add('active');\n", + " fig.buttons['Zoom'].classList.remove('active');\n", + " } else if (msg['mode'] === 'ZOOM') {\n", + " fig.buttons['Pan'].classList.remove('active');\n", + " fig.buttons['Zoom'].classList.add('active');\n", + " } else {\n", + " fig.buttons['Pan'].classList.remove('active');\n", + " fig.buttons['Zoom'].classList.remove('active');\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.updated_canvas_event = function () {\n", + " // Called whenever the canvas gets updated.\n", + " this.send_message('ack', {});\n", + "};\n", + "\n", + "// A function to construct a web socket function for onmessage handling.\n", + "// Called in the figure constructor.\n", + "mpl.figure.prototype._make_on_message_function = function (fig) {\n", + " return function socket_on_message(evt) {\n", + " if (evt.data instanceof Blob) {\n", + " var img = evt.data;\n", + " if (img.type !== 'image/png') {\n", + " /* FIXME: We get \"Resource interpreted as Image but\n", + " * transferred with MIME type text/plain:\" errors on\n", + " * Chrome. But how to set the MIME type? It doesn't seem\n", + " * to be part of the websocket stream */\n", + " img.type = 'image/png';\n", + " }\n", + "\n", + " /* Free the memory for the previous frames */\n", + " if (fig.imageObj.src) {\n", + " (window.URL || window.webkitURL).revokeObjectURL(\n", + " fig.imageObj.src\n", + " );\n", + " }\n", + "\n", + " fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n", + " img\n", + " );\n", + " fig.updated_canvas_event();\n", + " fig.waiting = false;\n", + " return;\n", + " } else if (\n", + " typeof evt.data === 'string' &&\n", + " evt.data.slice(0, 21) === 'data:image/png;base64'\n", + " ) {\n", + " fig.imageObj.src = evt.data;\n", + " fig.updated_canvas_event();\n", + " fig.waiting = false;\n", + " return;\n", + " }\n", + "\n", + " var msg = JSON.parse(evt.data);\n", + " var msg_type = msg['type'];\n", + "\n", + " // Call the \"handle_{type}\" callback, which takes\n", + " // the figure and JSON message as its only arguments.\n", + " try {\n", + " var callback = fig['handle_' + msg_type];\n", + " } catch (e) {\n", + " console.log(\n", + " \"No handler for the '\" + msg_type + \"' message type: \",\n", + " msg\n", + " );\n", + " return;\n", + " }\n", + "\n", + " if (callback) {\n", + " try {\n", + " // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n", + " callback(fig, msg);\n", + " } catch (e) {\n", + " console.log(\n", + " \"Exception inside the 'handler_\" + msg_type + \"' callback:\",\n", + " e,\n", + " e.stack,\n", + " msg\n", + " );\n", + " }\n", + " }\n", + " };\n", + "};\n", + "\n", + "// from https://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n", + "mpl.findpos = function (e) {\n", + " //this section is from http://www.quirksmode.org/js/events_properties.html\n", + " var targ;\n", + " if (!e) {\n", + " e = window.event;\n", + " }\n", + " if (e.target) {\n", + " targ = e.target;\n", + " } else if (e.srcElement) {\n", + " targ = e.srcElement;\n", + " }\n", + " if (targ.nodeType === 3) {\n", + " // defeat Safari bug\n", + " targ = targ.parentNode;\n", + " }\n", + "\n", + " // pageX,Y are the mouse positions relative to the document\n", + " var boundingRect = targ.getBoundingClientRect();\n", + " var x = e.pageX - (boundingRect.left + document.body.scrollLeft);\n", + " var y = e.pageY - (boundingRect.top + document.body.scrollTop);\n", + "\n", + " return { x: x, y: y };\n", + "};\n", + "\n", + "/*\n", + " * return a copy of an object with only non-object keys\n", + " * we need this to avoid circular references\n", + " * https://stackoverflow.com/a/24161582/3208463\n", + " */\n", + "function simpleKeys(original) {\n", + " return Object.keys(original).reduce(function (obj, key) {\n", + " if (typeof original[key] !== 'object') {\n", + " obj[key] = original[key];\n", + " }\n", + " return obj;\n", + " }, {});\n", + "}\n", + "\n", + "mpl.figure.prototype.mouse_event = function (event, name) {\n", + " var canvas_pos = mpl.findpos(event);\n", + "\n", + " if (name === 'button_press') {\n", + " this.canvas.focus();\n", + " this.canvas_div.focus();\n", + " }\n", + "\n", + " var x = canvas_pos.x * this.ratio;\n", + " var y = canvas_pos.y * this.ratio;\n", + "\n", + " this.send_message(name, {\n", + " x: x,\n", + " y: y,\n", + " button: event.button,\n", + " step: event.step,\n", + " guiEvent: simpleKeys(event),\n", + " });\n", + "\n", + " /* This prevents the web browser from automatically changing to\n", + " * the text insertion cursor when the button is pressed. We want\n", + " * to control all of the cursor setting manually through the\n", + " * 'cursor' event from matplotlib */\n", + " event.preventDefault();\n", + " return false;\n", + "};\n", + "\n", + "mpl.figure.prototype._key_event_extra = function (_event, _name) {\n", + " // Handle any extra behaviour associated with a key event\n", + "};\n", + "\n", + "mpl.figure.prototype.key_event = function (event, name) {\n", + " // Prevent repeat events\n", + " if (name === 'key_press') {\n", + " if (event.key === this._key) {\n", + " return;\n", + " } else {\n", + " this._key = event.key;\n", + " }\n", + " }\n", + " if (name === 'key_release') {\n", + " this._key = null;\n", + " }\n", + "\n", + " var value = '';\n", + " if (event.ctrlKey && event.key !== 'Control') {\n", + " value += 'ctrl+';\n", + " }\n", + " else if (event.altKey && event.key !== 'Alt') {\n", + " value += 'alt+';\n", + " }\n", + " else if (event.shiftKey && event.key !== 'Shift') {\n", + " value += 'shift+';\n", + " }\n", + "\n", + " value += 'k' + event.key;\n", + "\n", + " this._key_event_extra(event, name);\n", + "\n", + " this.send_message(name, { key: value, guiEvent: simpleKeys(event) });\n", + " return false;\n", + "};\n", + "\n", + "mpl.figure.prototype.toolbar_button_onclick = function (name) {\n", + " if (name === 'download') {\n", + " this.handle_save(this, null);\n", + " } else {\n", + " this.send_message('toolbar_button', { name: name });\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.toolbar_button_onmouseover = function (tooltip) {\n", + " this.message.textContent = tooltip;\n", + "};\n", + "\n", + "///////////////// REMAINING CONTENT GENERATED BY embed_js.py /////////////////\n", + "// prettier-ignore\n", + "var _JSXTOOLS_RESIZE_OBSERVER=function(A){var t,i=new WeakMap,n=new WeakMap,a=new WeakMap,r=new WeakMap,o=new Set;function s(e){if(!(this instanceof s))throw new TypeError(\"Constructor requires 'new' operator\");i.set(this,e)}function h(){throw new TypeError(\"Function is not a constructor\")}function c(e,t,i,n){e=0 in arguments?Number(arguments[0]):0,t=1 in arguments?Number(arguments[1]):0,i=2 in arguments?Number(arguments[2]):0,n=3 in arguments?Number(arguments[3]):0,this.right=(this.x=this.left=e)+(this.width=i),this.bottom=(this.y=this.top=t)+(this.height=n),Object.freeze(this)}function d(){t=requestAnimationFrame(d);var s=new WeakMap,p=new Set;o.forEach((function(t){r.get(t).forEach((function(i){var r=t instanceof window.SVGElement,o=a.get(t),d=r?0:parseFloat(o.paddingTop),f=r?0:parseFloat(o.paddingRight),l=r?0:parseFloat(o.paddingBottom),u=r?0:parseFloat(o.paddingLeft),g=r?0:parseFloat(o.borderTopWidth),m=r?0:parseFloat(o.borderRightWidth),w=r?0:parseFloat(o.borderBottomWidth),b=u+f,F=d+l,v=(r?0:parseFloat(o.borderLeftWidth))+m,W=g+w,y=r?0:t.offsetHeight-W-t.clientHeight,E=r?0:t.offsetWidth-v-t.clientWidth,R=b+v,z=F+W,M=r?t.width:parseFloat(o.width)-R-E,O=r?t.height:parseFloat(o.height)-z-y;if(n.has(t)){var k=n.get(t);if(k[0]===M&&k[1]===O)return}n.set(t,[M,O]);var S=Object.create(h.prototype);S.target=t,S.contentRect=new c(u,d,M,O),s.has(i)||(s.set(i,[]),p.add(i)),s.get(i).push(S)}))})),p.forEach((function(e){i.get(e).call(e,s.get(e),e)}))}return s.prototype.observe=function(i){if(i instanceof window.Element){r.has(i)||(r.set(i,new Set),o.add(i),a.set(i,window.getComputedStyle(i)));var n=r.get(i);n.has(this)||n.add(this),cancelAnimationFrame(t),t=requestAnimationFrame(d)}},s.prototype.unobserve=function(i){if(i instanceof window.Element&&r.has(i)){var n=r.get(i);n.has(this)&&(n.delete(this),n.size||(r.delete(i),o.delete(i))),n.size||r.delete(i),o.size||cancelAnimationFrame(t)}},A.DOMRectReadOnly=c,A.ResizeObserver=s,A.ResizeObserverEntry=h,A}; // eslint-disable-line\n", + "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Left button pans, Right button zooms\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-arrows\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\\nx/y fixes axis\", \"fa fa-square-o\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o\", \"download\"]];\n", + "\n", + "mpl.extensions = [\"eps\", \"jpeg\", \"pgf\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\", \"webp\"];\n", + "\n", + "mpl.default_extension = \"png\";/* global mpl */\n", + "\n", + "var comm_websocket_adapter = function (comm) {\n", + " // Create a \"websocket\"-like object which calls the given IPython comm\n", + " // object with the appropriate methods. Currently this is a non binary\n", + " // socket, so there is still some room for performance tuning.\n", + " var ws = {};\n", + "\n", + " ws.binaryType = comm.kernel.ws.binaryType;\n", + " ws.readyState = comm.kernel.ws.readyState;\n", + " function updateReadyState(_event) {\n", + " if (comm.kernel.ws) {\n", + " ws.readyState = comm.kernel.ws.readyState;\n", + " } else {\n", + " ws.readyState = 3; // Closed state.\n", + " }\n", + " }\n", + " comm.kernel.ws.addEventListener('open', updateReadyState);\n", + " comm.kernel.ws.addEventListener('close', updateReadyState);\n", + " comm.kernel.ws.addEventListener('error', updateReadyState);\n", + "\n", + " ws.close = function () {\n", + " comm.close();\n", + " };\n", + " ws.send = function (m) {\n", + " //console.log('sending', m);\n", + " comm.send(m);\n", + " };\n", + " // Register the callback with on_msg.\n", + " comm.on_msg(function (msg) {\n", + " //console.log('receiving', msg['content']['data'], msg);\n", + " var data = msg['content']['data'];\n", + " if (data['blob'] !== undefined) {\n", + " data = {\n", + " data: new Blob(msg['buffers'], { type: data['blob'] }),\n", + " };\n", + " }\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", + " ws.onmessage(data);\n", + " });\n", + " return ws;\n", + "};\n", + "\n", + "mpl.mpl_figure_comm = function (comm, msg) {\n", + " // This is the function which gets called when the mpl process\n", + " // starts-up an IPython Comm through the \"matplotlib\" channel.\n", + "\n", + " var id = msg.content.data.id;\n", + " // Get hold of the div created by the display call when the Comm\n", + " // socket was opened in Python.\n", + " var element = document.getElementById(id);\n", + " var ws_proxy = comm_websocket_adapter(comm);\n", + "\n", + " function ondownload(figure, _format) {\n", + " window.open(figure.canvas.toDataURL());\n", + " }\n", + "\n", + " var fig = new mpl.figure(id, ws_proxy, ondownload, element);\n", + "\n", + " // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n", + " // web socket which is closed, not our websocket->open comm proxy.\n", + " ws_proxy.onopen();\n", + "\n", + " fig.parent_element = element;\n", + " fig.cell_info = mpl.find_output_cell(\"
\");\n", + " if (!fig.cell_info) {\n", + " console.error('Failed to find cell for figure', id, fig);\n", + " return;\n", + " }\n", + " fig.cell_info[0].output_area.element.on(\n", + " 'cleared',\n", + " { fig: fig },\n", + " fig._remove_fig_handler\n", + " );\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_close = function (fig, msg) {\n", + " var width = fig.canvas.width / fig.ratio;\n", + " fig.cell_info[0].output_area.element.off(\n", + " 'cleared',\n", + " fig._remove_fig_handler\n", + " );\n", + " fig.resizeObserverInstance.unobserve(fig.canvas_div);\n", + "\n", + " // Update the output cell to use the data from the current canvas.\n", + " fig.push_to_output();\n", + " var dataURL = fig.canvas.toDataURL();\n", + " // Re-enable the keyboard manager in IPython - without this line, in FF,\n", + " // the notebook keyboard shortcuts fail.\n", + " IPython.keyboard_manager.enable();\n", + " fig.parent_element.innerHTML =\n", + " '';\n", + " fig.close_ws(fig, msg);\n", + "};\n", + "\n", + "mpl.figure.prototype.close_ws = function (fig, msg) {\n", + " fig.send_message('closing', msg);\n", + " // fig.ws.close()\n", + "};\n", + "\n", + "mpl.figure.prototype.push_to_output = function (_remove_interactive) {\n", + " // Turn the data on the canvas into data in the output cell.\n", + " var width = this.canvas.width / this.ratio;\n", + " var dataURL = this.canvas.toDataURL();\n", + " this.cell_info[1]['text/html'] =\n", + " '';\n", + "};\n", + "\n", + "mpl.figure.prototype.updated_canvas_event = function () {\n", + " // Tell IPython that the notebook contents must change.\n", + " IPython.notebook.set_dirty(true);\n", + " this.send_message('ack', {});\n", + " var fig = this;\n", + " // Wait a second, then push the new image to the DOM so\n", + " // that it is saved nicely (might be nice to debounce this).\n", + " setTimeout(function () {\n", + " fig.push_to_output();\n", + " }, 1000);\n", + "};\n", + "\n", + "mpl.figure.prototype._init_toolbar = function () {\n", + " var fig = this;\n", + "\n", + " var toolbar = document.createElement('div');\n", + " toolbar.classList = 'btn-toolbar';\n", + " this.root.appendChild(toolbar);\n", + "\n", + " function on_click_closure(name) {\n", + " return function (_event) {\n", + " return fig.toolbar_button_onclick(name);\n", + " };\n", + " }\n", + "\n", + " function on_mouseover_closure(tooltip) {\n", + " return function (event) {\n", + " if (!event.currentTarget.disabled) {\n", + " return fig.toolbar_button_onmouseover(tooltip);\n", + " }\n", + " };\n", + " }\n", + "\n", + " fig.buttons = {};\n", + " var buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'btn-group';\n", + " var button;\n", + " for (var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " /* Instead of a spacer, we start a new button group. */\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + " buttonGroup = document.createElement('div');\n", + " buttonGroup.classList = 'btn-group';\n", + " continue;\n", + " }\n", + "\n", + " button = fig.buttons[name] = document.createElement('button');\n", + " button.classList = 'btn btn-default';\n", + " button.href = '#';\n", + " button.title = name;\n", + " button.innerHTML = '';\n", + " button.addEventListener('click', on_click_closure(method_name));\n", + " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", + " buttonGroup.appendChild(button);\n", + " }\n", + "\n", + " if (buttonGroup.hasChildNodes()) {\n", + " toolbar.appendChild(buttonGroup);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = document.createElement('span');\n", + " status_bar.classList = 'mpl-message pull-right';\n", + " toolbar.appendChild(status_bar);\n", + " this.message = status_bar;\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = document.createElement('div');\n", + " buttongrp.classList = 'btn-group inline pull-right';\n", + " button = document.createElement('button');\n", + " button.classList = 'btn btn-mini btn-primary';\n", + " button.href = '#';\n", + " button.title = 'Stop Interaction';\n", + " button.innerHTML = '';\n", + " button.addEventListener('click', function (_evt) {\n", + " fig.handle_close(fig, {});\n", + " });\n", + " button.addEventListener(\n", + " 'mouseover',\n", + " on_mouseover_closure('Stop Interaction')\n", + " );\n", + " buttongrp.appendChild(button);\n", + " var titlebar = this.root.querySelector('.ui-dialog-titlebar');\n", + " titlebar.insertBefore(buttongrp, titlebar.firstChild);\n", + "};\n", + "\n", + "mpl.figure.prototype._remove_fig_handler = function (event) {\n", + " var fig = event.data.fig;\n", + " if (event.target !== this) {\n", + " // Ignore bubbled events from children.\n", + " return;\n", + " }\n", + " fig.close_ws(fig, {});\n", + "};\n", + "\n", + "mpl.figure.prototype._root_extra_style = function (el) {\n", + " el.style.boxSizing = 'content-box'; // override notebook setting of border-box.\n", + "};\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function (el) {\n", + " // this is important to make the div 'focusable\n", + " el.setAttribute('tabindex', 0);\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " } else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype._key_event_extra = function (event, _name) {\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which === 13) {\n", + " this.canvas_div.blur();\n", + " // select the cell after this one\n", + " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", + " IPython.notebook.select(index + 1);\n", + " }\n", + "};\n", + "\n", + "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", + " fig.ondownload(fig, null);\n", + "};\n", + "\n", + "mpl.find_output_cell = function (html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i = 0; i < ncells; i++) {\n", + " var cell = cells[i];\n", + " if (cell.cell_type === 'code') {\n", + " for (var j = 0; j < cell.output_area.outputs.length; j++) {\n", + " var data = cell.output_area.outputs[j];\n", + " if (data.data) {\n", + " // IPython >= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] === html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "};\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel !== null) {\n", + " IPython.notebook.kernel.comm_manager.register_target(\n", + " 'matplotlib',\n", + " mpl.mpl_figure_comm\n", + " );\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(ncols=3, figsize=(9, 3), sharey=True, sharex=True)\n", + "fig.subplots_adjust(wspace=0)\n", + "for m, k in enumerate(ids):\n", + " r200 = cat[n][\"r200\"][k]\n", + " m200 = cat[n][\"m200\"][k]\n", + " \n", + " axs[m].set_title(r\"$\\log M_{{200c}} / M_\\odot = {:.3f}$\".format(np.log10(m200)))\n", + " for i, j in enumerate(matcher.search_sim_indices(n)):\n", + " indxs, dist = match[i]\n", + " indxs = indxs[k]\n", + " dist = dist[k]\n", + "\n", + " dlogmass = np.abs(np.log10(cat[j][\"m200\"][indxs] / m200))\n", + " axs[m].scatter(dist / r200, dlogmass)\n", + " \n", + " \n", + "for i in range(3):\n", + " axs[i].set_xlabel(r\"$\\Delta r_i / R_{200c}$\")\n", + "axs[0].set_ylabel(r\"$|\\log \\dfrac{M_i}{M_{200c}}|$\")\n", + "plt.tight_layout(w_pad=0)\n", + "# plt.savefig(\"../plots/lowest_massive.png\", dpi=450)\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f0f59246", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a25c7146", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2930abd", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T10:56:11.296860Z", + "start_time": "2022-11-21T10:56:10.813212Z" + } + }, + "outputs": [], + "source": [ + "plt.figure()\n", + "for k in range(cat.N-1):\n", + " plt.scatter(x[k, :, 1], x[k, :, 2])\n", + "\n", + " \n", + "plt.xlabel(r\"$\\Delta r_i / R_{200c}$\")\n", + "plt.ylabel(r\"$|\\log \\dfrac{M_i}{M_{200c}}|$\")\n", + "plt.savefig(\"../plots/example.png\", dpi=450)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5101f569", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T10:36:59.128848Z", + "start_time": "2022-11-21T10:36:58.543085Z" + } + }, + "outputs": [], + "source": [ + "x[:, 0, 2]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6794c154", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa252ad3", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T07:38:33.029175Z", + "start_time": "2022-11-21T07:38:32.649968Z" + } + }, + "outputs": [], + "source": [ + "cat.cross_knn_position(0, 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8dc14292", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f7cd638", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d655f5de", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9894fd1a", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41c21c3c", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T07:05:14.278302Z", + "start_time": "2022-11-21T07:05:14.249214Z" + } + }, + "outputs": [], + "source": [ + "x = cat.positions[6, :].reshape(-1, 3)\n", + "\n", + "\n", + "# x = np.asarray([0., 0., 0.]).reshape(-1, 3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "440c00a3", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T07:07:55.025995Z", + "start_time": "2022-11-21T07:07:54.988489Z" + } + }, + "outputs": [], + "source": [ + "cat.knn_position(x, 5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c510eca", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c59341c0", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:55:17.697688Z", + "start_time": "2022-11-21T06:55:17.667672Z" + } + }, + "outputs": [], + "source": [ + "dist, knns = cat.knn_position([5, 9], 5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "004c5a79", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:55:40.505297Z", + "start_time": "2022-11-21T06:55:40.475623Z" + } + }, + "outputs": [], + "source": [ + "dist" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "efcc87ac", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:55:47.645798Z", + "start_time": "2022-11-21T06:55:47.615986Z" + } + }, + "outputs": [], + "source": [ + "knns" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bf75953", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:27:19.883364Z", + "start_time": "2022-11-21T06:27:19.435813Z" + } + }, + "outputs": [], + "source": [ + "model = NearestNeighbors()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e751553", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:27:38.874030Z", + "start_time": "2022-11-21T06:27:38.663757Z" + } + }, + "outputs": [], + "source": [ + "model.fit(cat[0].positions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8cb2bcf", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:37:19.127300Z", + "start_time": "2022-11-21T06:37:18.905607Z" + } + }, + "outputs": [], + "source": [ + "p = cat[0].positions[:2, :].reshape(-1, 3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae96df81", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:37:19.354558Z", + "start_time": "2022-11-21T06:37:19.324845Z" + } + }, + "outputs": [], + "source": [ + "p" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e7b3655", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-21T06:37:20.172812Z", + "start_time": "2022-11-21T06:37:20.141347Z" + } + }, + "outputs": [], + "source": [ + "model.kneighbors(p, n_neighbors=2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d3bcb85", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "307f3b62", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a885af97", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T16:47:33.606111Z", + "start_time": "2022-11-20T16:47:29.246847Z" + } + }, + "outputs": [], + "source": [ + "cat = csiborgtools.io.HaloCatalogue(9844, 1016, minimum_m500=0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b73c9bc", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T16:47:52.218696Z", + "start_time": "2022-11-20T16:47:52.097392Z" + } + }, + "outputs": [], + "source": [ + "cat[\"dec\"].size" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd76ed7a", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22abdd5c", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T15:15:05.423239Z", + "start_time": "2022-11-20T15:15:05.366755Z" + } + }, + "outputs": [], + "source": [ + "planck = csiborgtools.io.PlanckClusters(\"../data/HFI_PCCS_SZ-union_R2.08.fits\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20855f15", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T15:15:05.479763Z", + "start_time": "2022-11-20T15:15:05.425076Z" + } + }, + "outputs": [], + "source": [ + "mcxc = csiborgtools.io.MCXCClusters(\"../data/mcxc.fits\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60e7d566", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T15:15:38.743815Z", + "start_time": "2022-11-20T15:15:38.663202Z" + } + }, + "outputs": [], + "source": [ + "planck.match_to_mcxc(mcxc)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2941d725", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T14:36:05.559176Z", + "start_time": "2022-11-20T14:36:05.143516Z" + } + }, + "outputs": [], + "source": [ + "planck[\"MSZ\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6347d993", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43293c40", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T14:34:56.350019Z", + "start_time": "2022-11-20T14:34:56.292284Z" + } + }, + "outputs": [], + "source": [ + "planck.data.dtype.names" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "509d4fbb", + "metadata": { + "ExecuteTime": { + "end_time": "2022-11-20T17:00:49.100716Z", + "start_time": "2022-11-20T17:00:48.529803Z" + } + }, + "outputs": [], + "source": [ + "csiborgtools.io.get_csiborg_ids(\"/mnt/extraspace/hdesmond\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, "id": "f39874e4", "metadata": { "ExecuteTime": { @@ -47,7 +2647,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "d243cc59", "metadata": { "ExecuteTime": { @@ -65,7 +2665,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "344e6b96", "metadata": { "ExecuteTime": { @@ -80,7 +2680,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "id": "d8bb7dce", "metadata": { "ExecuteTime": { @@ -95,7 +2695,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, "id": "71026a85", "metadata": { "ExecuteTime": { @@ -110,7 +2710,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "id": "3ca4a7b7", "metadata": { "ExecuteTime": { @@ -118,28 +2718,14 @@ "start_time": "2022-11-20T12:51:52.659741Z" } }, - "outputs": [ - { - "data": { - "text/html": [ - "
NearestNeighbors()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
" - ], - "text/plain": [ - "NearestNeighbors()" - ] - }, - "execution_count": 14, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "neighbors.fit(X)" ] }, { "cell_type": "code", - "execution_count": 19, + "execution_count": null, "id": "6e5ec9c8", "metadata": { "ExecuteTime": { @@ -147,19 +2733,7 @@ "start_time": "2022-11-20T12:52:50.429353Z" } }, - "outputs": [ - { - "data": { - "text/plain": [ - "(array([[0. , 0.00441153, 0.00459373, 0.00568435, 0.00596298]]),\n", - " array([[ 0, 186, 513, 130, 1795]]))" - ] - }, - "execution_count": 19, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "p = X[0, :]\n", "\n", @@ -168,7 +2742,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "id": "304d53aa", "metadata": { "ExecuteTime": { @@ -176,23 +2750,12 @@ "start_time": "2022-11-20T12:52:31.796380Z" } }, - "outputs": [ - { - "data": { - "text/plain": [ - "array([-0.04183578, -0.21976852, -0.05943659])" - ] - }, - "execution_count": 16, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "a65dac24", "metadata": { "ExecuteTime": { @@ -200,19 +2763,7 @@ "start_time": "2022-11-20T12:49:09.416602Z" } }, - "outputs": [ - { - "data": { - "text/plain": [ - "array([-0.21976852, -0.15894653, -0.16384361, ..., -0.23763191,\n", - " -0.23351993, -0.23352051])" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "data[\"peak_y\"]" ] diff --git a/scripts/run_fit_halos.py b/scripts/run_fit_halos.py index 6bb084d..7c8008e 100644 --- a/scripts/run_fit_halos.py +++ b/scripts/run_fit_halos.py @@ -32,9 +32,6 @@ import utils F64 = numpy.float64 I64 = numpy.int64 -# Simulations and their snapshot to analyze -Nsims = [9844] -Nsnap = 1016 # Get MPI things comm = MPI.COMM_WORLD @@ -50,68 +47,76 @@ cols_collect = [("npart", I64), ("totpartmass", F64), ("Rs", F64), ("rmax", F64), ("r200", F64), ("r500", F64), ("m200", F64), ("m500", F64)] -# NOTE later loop over sims too -Nsim = Nsims[0] -simpath = csiborgtools.io.get_sim_path(Nsim) -box = csiborgtools.units.BoxUnits(Nsnap, simpath) +Nsims = csiborgtools.read.get_csiborg_ids("/mnt/extraspace/hdesmond") +for i, Nsim in enumerate(Nsims): + if rank == 0: + print("{}: calculating {}th simulation.".format(datetime.now(), i)) -jobs = csiborgtools.fits.split_jobs(utils.Nsplits, nproc)[rank] -for icount, Nsplit in enumerate(jobs): - print("{}: rank {} working {} / {} jobs.".format(datetime.now(), rank, - icount + 1, len(jobs))) - parts, part_clumps, clumps = csiborgtools.fits.load_split_particles( - Nsplit, loaddir, Nsim, Nsnap, remove_split=False) + simpath = csiborgtools.read.get_sim_path(Nsim) + Nsnap = csiborgtools.read.get_maximum_snapshot(simpath) + box = csiborgtools.units.BoxUnits(Nsnap, simpath) - N = clumps.size - cols = [("index", I64), ("npart", I64), ("totpartmass", F64), - ("Rs", F64), ("rho0", F64), ("conc", F64), - ("vx", F64), ("vy", F64), ("vz", F64), - ("rmin", F64), ("rmax", F64), - ("r200", F64), ("r500", F64), ("m200", F64), ("m500", F64)] - out = csiborgtools.utils.cols_to_structured(N, cols) - out["index"] = clumps["index"] + jobs = csiborgtools.fits.split_jobs(utils.Nsplits, nproc)[rank] + for Nsplit in jobs: + parts, part_clumps, clumps = csiborgtools.fits.load_split_particles( + Nsplit, loaddir, Nsim, Nsnap, remove_split=False) - for n in range(N): - # Pick clump and its particles - xs = csiborgtools.fits.pick_single_clump(n, parts, part_clumps, clumps) - clump = csiborgtools.fits.Clump.from_arrays(*xs, rhoc=box.box_rhoc) - out["npart"][n] = clump.Npart - out["rmin"][n] = clump.rmin - out["rmax"][n] = clump.rmax - out["totpartmass"][n] = clump.total_particle_mass - out["vx"] = numpy.average(clump.vel[:, 0], weights=clump.m) - out["vy"] = numpy.average(clump.vel[:, 1], weights=clump.m) - out["vz"] = numpy.average(clump.vel[:, 2], weights=clump.m) + N = clumps.size + cols = [("index", I64), ("npart", I64), ("totpartmass", F64), + ("Rs", F64), ("rho0", F64), ("conc", F64), + ("vx", F64), ("vy", F64), ("vz", F64), + ("rmin", F64), ("rmax", F64), + ("r200", F64), ("r500", F64), ("m200", F64), ("m500", F64)] + out = csiborgtools.utils.cols_to_structured(N, cols) + out["index"] = clumps["index"] - # Spherical overdensity radii and masses - rs, ms = clump.spherical_overdensity_mass([200, 500]) - out["r200"][n] = rs[0] - out["r500"][n] = rs[1] - out["m200"][n] = ms[0] - out["m500"][n] = ms[1] + for n in range(N): + # Pick clump and its particles + xs = csiborgtools.fits.pick_single_clump(n, parts, part_clumps, + clumps) + clump = csiborgtools.fits.Clump.from_arrays(*xs, rhoc=box.box_rhoc) + out["npart"][n] = clump.Npart + out["rmin"][n] = clump.rmin + out["rmax"][n] = clump.rmax + out["totpartmass"][n] = clump.total_particle_mass + out["vx"] = numpy.average(clump.vel[:, 0], weights=clump.m) + out["vy"] = numpy.average(clump.vel[:, 1], weights=clump.m) + out["vz"] = numpy.average(clump.vel[:, 2], weights=clump.m) - # NFW profile fit - if clump.Npart > 10 and numpy.isfinite(out["r200"][n]): - nfwpost = csiborgtools.fits.NFWPosterior(clump) - logRs, __ = nfwpost.maxpost_logRs() - Rs = 10**logRs - if not numpy.isnan(logRs): - out["Rs"][n] = Rs - out["rho0"][n] = nfwpost.rho0_from_Rs(Rs) - out["conc"][n] = out["r200"][n] / Rs + # Spherical overdensity radii and masses + rs, ms = clump.spherical_overdensity_mass([200, 500]) + out["r200"][n] = rs[0] + out["r500"][n] = rs[1] + out["m200"][n] = ms[0] + out["m500"][n] = ms[1] - csiborgtools.io.dump_split(out, Nsplit, Nsim, Nsnap, dumpdir) + # NFW profile fit + if clump.Npart > 10 and numpy.isfinite(out["r200"][n]): + nfwpost = csiborgtools.fits.NFWPosterior(clump) + logRs, __ = nfwpost.maxpost_logRs() + Rs = 10**logRs + if not numpy.isnan(logRs): + out["Rs"][n] = Rs + out["rho0"][n] = nfwpost.rho0_from_Rs(Rs) + out["conc"][n] = out["r200"][n] / Rs + + csiborgtools.read.dump_split(out, Nsplit, Nsim, Nsnap, dumpdir) + + # Wait until all jobs finished before moving to another simulation + comm.Barrier() + + # Use the rank 0 to combine outputs for this CSiBORG realisation + if rank == 0: + print("Collecting results!") + out_collected = csiborgtools.read.combine_splits( + utils.Nsplits, Nsim, Nsnap, utils.dumpdir, cols_collect, + remove_splits=True, verbose=False) + fname = join(utils.dumpdir, "ramses_out_{}_{}.npy" + .format(str(Nsim).zfill(5), str(Nsnap).zfill(5))) + print("Saving results to `{}`.".format(fname)) + numpy.save(fname, out_collected) + + comm.Barrier() -# Force all ranks to wait -comm.Barrier() -# Use the rank 0 to combine outputs for this CSiBORG realisation if rank == 0: - print("Collecting results!") - out_collected = csiborgtools.io.combine_splits( - utils.Nsplits, Nsim, Nsnap, utils.dumpdir, cols_collect, - remove_splits=True, verbose=False) - fname = join(utils.dumpdir, "ramses_out_{}_{}.npy" - .format(str(Nsim).zfill(5), str(Nsnap).zfill(5))) - print("Saving results to `{}`.".format(fname)) - numpy.save(fname, out_collected) print("All finished! See ya!") diff --git a/scripts/run_split_halos.py b/scripts/run_split_halos.py index fe7589b..ae33581 100644 --- a/scripts/run_split_halos.py +++ b/scripts/run_split_halos.py @@ -34,7 +34,7 @@ comm = MPI.COMM_WORLD rank = comm.Get_rank() nproc = comm.Get_size() -Nsims = csiborgtools.io.get_csiborg_ids("/mnt/extraspace/hdesmond") +Nsims = csiborgtools.read.get_csiborg_ids("/mnt/extraspace/hdesmond") partcols = ["x", "y", "z", "vx", "vy", "vz", "M", "level"] dumpdir = join(utils.dumpdir, "temp") @@ -43,16 +43,16 @@ for icount, sim_index in enumerate(jobs): print("{}: rank {} working {} / {} jobs.".format(datetime.now(), rank, icount + 1, len(jobs))) Nsim = Nsims[sim_index] - simpath = csiborgtools.io.get_sim_path(Nsim) - Nsnap = csiborgtools.io.get_maximum_snapshot(simpath) + simpath = csiborgtools.read.get_sim_path(Nsim) + Nsnap = csiborgtools.read.get_maximum_snapshot(simpath) # Load the clumps, particles' clump IDs and particles. - clumps = csiborgtools.io.read_clumps(Nsnap, simpath) - particle_clumps = csiborgtools.io.read_clumpid(Nsnap, simpath, - verbose=False) - particles = csiborgtools.io.read_particle(partcols, Nsnap, simpath, - verbose=False) + clumps = csiborgtools.read.read_clumps(Nsnap, simpath) + particle_clumps = csiborgtools.read.read_clumpid(Nsnap, simpath, + verbose=False) + particles = csiborgtools.read.read_particle(partcols, Nsnap, simpath, + verbose=False) # Drop all particles whose clump index is 0 (not assigned to any halo) - particle_clumps, particles = csiborgtools.io.drop_zero_indx( + particle_clumps, particles = csiborgtools.read.drop_zero_indx( particle_clumps, particles) # Dump it! csiborgtools.fits.dump_split_particles(particles, particle_clumps, clumps, diff --git a/scripts/utils.py b/scripts/utils.py index 5d88c3c..4b64e88 100644 --- a/scripts/utils.py +++ b/scripts/utils.py @@ -16,16 +16,14 @@ Notebook utility functions. """ +# import numpy +# from os.path import join -import numpy -from os.path import join -from astropy.cosmology import FlatLambdaCDM - -try: - import csiborgtools -except ModuleNotFoundError: - import sys - sys.path.append("../") +# try: +# import csiborgtools +# except ModuleNotFoundError: +# import sys +# sys.path.append("../") Nsplits = 200 @@ -42,52 +40,3 @@ _virgo = {"RA": (12 + 27 / 60) * 15, "COMDIST": 16.5} specific_clusters = {"Coma": _coma, "Virgo": _virgo} - - -def load_processed(Nsim, Nsnap): - simpath = csiborgtools.io.get_sim_path(Nsim) - outfname = join( - dumpdir, "ramses_out_{}_{}.npy".format(str(Nsim).zfill(5), - str(Nsnap).zfill(5))) - data = numpy.load(outfname) - # Add mmain - mmain = csiborgtools.io.read_mmain(Nsim, "/mnt/zfsusers/hdesmond/Mmain") - data = csiborgtools.io.merge_mmain_to_clumps(data, mmain) - csiborgtools.utils.flip_cols(data, "peak_x", "peak_z") - # Cut on numbre of particles and finite m200 - data = data[(data["npart"] > 100) & numpy.isfinite(data["m200"])] - - # Do unit conversion - boxunits = csiborgtools.units.BoxUnits(Nsnap, simpath) - convert_cols = ["m200", "m500", "totpartmass", "mass_mmain", - "r200", "r500", "Rs", "rho0", "peak_x", "peak_y", "peak_z"] - data = csiborgtools.units.convert_from_boxunits( - data, convert_cols, boxunits) - # Now calculate spherical coordinates - d, ra, dec = csiborgtools.units.cartesian_to_radec(data) - data = csiborgtools.utils.add_columns( - data, [d, ra, dec], ["dist", "ra", "dec"]) - return data, boxunits - - -def load_planck2015(max_comdist=214): - cosmo = FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728) - fpath = "../data/HFI_PCCS_SZ-union_R2.08.fits" - return csiborgtools.io.read_planck2015(fpath, cosmo, max_comdist) - - -def load_mcxc(max_comdist=214): - cosmo = FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728) - fpath = ("../data/mcxc.fits") - return csiborgtools.io.read_mcxc(fpath, cosmo, max_comdist) - - -def load_2mpp(): - cosmo = FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728) - return csiborgtools.io.read_2mpp("../data/2M++_galaxy_catalog.dat", cosmo) - - -def load_2mpp_groups(): - cosmo = FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728) - return csiborgtools.io.read_2mpp_groups( - "../data/../data/2M++_group_catalog.dat", cosmo)