diff --git a/csiborgtools/read/__init__.py b/csiborgtools/read/__init__.py index 9e16e94..026275a 100644 --- a/csiborgtools/read/__init__.py +++ b/csiborgtools/read/__init__.py @@ -21,9 +21,9 @@ from .obs import ( # noqa TwoMPPGalaxies, TwoMPPGroups, ) -from .outsim import combine_splits, dump_split # noqa from .overlap_summary import NPairsOverlap, PairOverlap, binned_resample_mean # noqa from .paths import CSiBORGPaths # noqa from .pk_summary import PKReader # noqa from .readsim import MmainReader, ParticleReader, halfwidth_select, read_initcm # noqa from .tpcf_summary import TPCFReader # noqa +from .utils import cartesian_to_radec, radec_to_cartesian \ No newline at end of file diff --git a/csiborgtools/read/halo_cat.py b/csiborgtools/read/halo_cat.py index af2d7e5..149648c 100644 --- a/csiborgtools/read/halo_cat.py +++ b/csiborgtools/read/halo_cat.py @@ -20,8 +20,8 @@ from sklearn.neighbors import NearestNeighbors from .box_units import BoxUnits from .paths import CSiBORGPaths -from .readsim import ParticleReader, read_initcm -from .utils import add_columns, cartesian_to_radec, flip_cols +from .readsim import ParticleReader +from .utils import cartesian_to_radec, flip_cols, radec_to_cartesian class BaseCatalogue(ABC): @@ -47,7 +47,6 @@ class BaseCatalogue(ABC): @nsim.setter def nsim(self, nsim): - assert isinstance(nsim, int) self._nsim = nsim @property @@ -136,11 +135,10 @@ class BaseCatalogue(ABC): ps = ['x0', 'y0', 'z0'] else: ps = ['x', 'y', 'z'] - pos = [self[p] for p in ps] - if cartesian: - return numpy.vstack(pos).T - else: - return numpy.vstack([cartesian_to_radec(*pos)]).T + pos = numpy.vstack([self[p] for p in ps]).T + if not cartesian: + pos = cartesian_to_radec(pos) + return pos def velocity(self): """ @@ -207,6 +205,64 @@ class BaseCatalogue(ABC): knn = self.knn(in_initial) return knn.radius_neighbors(X, radius, sort_results=True) + def angular_neighbours(self, X, ang_radius, rad_tolerance=None): + r""" + Find nearest neighbours within `ang_radius` of query points `X`. + Optionally applies radial tolerance, which is expected to be in + :math:`\mathrm{cMpc}`. + + Parameters + ---------- + X : 2-dimensional array of shape `(n_queries, 2)` or `(n_queries, 3)` + Query positions. If 2-dimensional, then RA and DEC in degrees. + If 3-dimensional, then radial distance in :math:`\mathrm{cMpc}`, + RA and DEC in degrees. + ang_radius : float + Angular radius in degrees. + rad_tolerance : float, optional + Radial tolerance in :math:`\mathrm{cMpc}`. + + Returns + ------- + dist : array of 1-dimensional arrays of shape `(n_neighbours,)` + Distance of each neighbour to the query point. + ind : array of 1-dimensional arrays of shape `(n_neighbours,)` + Indices of each neighbour in this catalogue. + """ + assert X.ndim == 2 + # We first get positions of haloes in this catalogue, store their + # radial distance and normalise them to unit vectors. + pos = self.position(in_initial=False, cartesian=True) + raddist = numpy.linalg.norm(pos, axis=1) + pos /= raddist.reshape(-1, 1) + # We convert RAdec query positions to unit vectors. If no radial + # distance provided add it. + if X.shape[1] == 2: + X = numpy.vstack([numpy.ones_like(X[:, 0]), X[:, 0], X[:, 1]]).T + radquery = None + else: + radquery = X[:, 0] + + X = radec_to_cartesian(X) + knn = NearestNeighbors(metric="cosine") + knn.fit(pos) + # Convert angular radius to cosine difference. + metric_maxdist = 1 - numpy.cos(numpy.deg2rad(ang_radius)) + dist, ind = knn.radius_neighbors(X, radius=metric_maxdist, + sort_results=True) + # And the cosine difference to angular distance. + for i in range(X.shape[0]): + dist[i] = numpy.rad2deg(numpy.arccos(1 - dist[i])) + + # Apply the radial tolerance + if rad_tolerance is not None: + assert radquery is not None + for i in range(X.shape[0]): + mask = numpy.abs(raddist[ind[i]] - radquery) < rad_tolerance + dist[i] = dist[i][mask] + ind[i] = ind[i][mask] + return dist, ind + @property def keys(self): """Catalogue keys.""" @@ -240,8 +296,12 @@ class ClumpsCatalogue(BaseCatalogue): The maximum comoving distance of a halo. By default :math:`155.5 / 0.705 ~ \mathrm{Mpc}` with assumed :math:`h = 0.705`, which corresponds to the high-resolution region. + minmass : len-2 tuple, optional + Minimum mass. The first element is the catalogue key and the second is + the value. """ - def __init__(self, nsim, paths, maxdist=155.5 / 0.705): + def __init__(self, nsim, paths, maxdist=155.5 / 0.705, + minmass=("mass_cl", 1e12)): self.nsim = nsim self.paths = paths # Read in the clumps from the final snapshot @@ -259,6 +319,7 @@ class ClumpsCatalogue(BaseCatalogue): data = self.box.convert_from_boxunits(data, ['x', 'y', 'z', "mass_cl"]) mask = numpy.sqrt(data['x']**2 + data['y']**2 + data['z']**2) < maxdist + mask &= data[minmass[0]] > minmass[1] self._data = data[mask] @property @@ -272,117 +333,6 @@ class ClumpsCatalogue(BaseCatalogue): """ return self["index"] == self["parent"] - def _set_data(self, min_mass, max_dist, load_init): - """ - TODO: old later remove. - Loads the data, merges with mmain, does various coordinate transforms. - """ - # Load the processed data - data = numpy.load(self.paths.hcat_path(self.nsim)) - - # Load the mmain file and add it to the data - # TODO: read the mmain here -# mmain = read_mmain(self.nsim, self.paths.mmain_dir) -# data = self.merge_mmain_to_clumps(data, mmain) - flip_cols(data, "peak_x", "peak_z") - - # Cut on number of particles and finite m200. Do not change! Hardcoded - data = data[(data["npart"] > 100) & numpy.isfinite(data["m200"])] - - # Now also load the initial positions - if load_init: - - initcm = read_initcm(self.nsim, - self.paths.initmatch_path(self.nsim, "cm")) - if initcm is not None: - data = self.merge_initmatch_to_clumps(data, initcm) - flip_cols(data, "x0", "z0") - - # 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) - - # And do the unit transform - if load_init and initcm is not None: - data = self.box.convert_from_boxunits( - data, ["x0", "y0", "z0", "lagpatch"]) - - # Convert all that is not an integer to float32 - names = list(data.dtype.names) - formats = [] - for name in names: - if data[name].dtype.char in numpy.typecodes["AllInteger"]: - formats.append(numpy.int32) - else: - formats.append(numpy.float32) - dtype = numpy.dtype({"names": names, "formats": formats}) - - # Apply cuts on distance and total particle mass if any - data = data[data["dist"] < max_dist] if max_dist is not None else data - data = (data[data["totpartmass"] > min_mass] - if min_mass is not None else data) - - self._data = data.astype(dtype) - - def merge_mmain_to_clumps(self, clumps, mmain): - """ - TODO: old, later remove. - 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"]) - - def merge_initmatch_to_clumps(self, clumps, initcat): - """ - TODO: old, later remove. - Merge columns from the `init_cm` files to the `clump` file. - - Parameters - ---------- - clumps : structured array - Clumps structured array. - initcat : structured array - Catalog with the clumps initial centre of mass at z = 70. - - Returns - ------- - out : structured array - """ - # There are more initcat clumps, so check which ones have z = 0 - # and then downsample - mask = numpy.isin(initcat["ID"], clumps["index"]) - initcat = initcat[mask] - # Now the index ordering should match - if not numpy.alltrue(initcat["ID"] == clumps["index"]): - raise ValueError( - "Ordering of `initcat` and `clumps` is inconsistent.") - - X = numpy.full((clumps.size, 4), numpy.nan) - for i, p in enumerate(['x', 'y', 'z', "lagpatch"]): - X[:, i] = initcat[p] - return add_columns(clumps, X, ["x0", "y0", "z0", "lagpatch"]) - class HaloCatalogue(BaseCatalogue): r""" @@ -403,8 +353,12 @@ class HaloCatalogue(BaseCatalogue): The maximum comoving distance of a halo. By default :math:`155.5 / 0.705 ~ \mathrm{Mpc}` with assumed :math:`h = 0.705`, which corresponds to the high-resolution region. + minmass : len-2 tuple + Minimum mass. The first element is the catalogue key and the second is + the value. """ - def __init__(self, nsim, paths, maxdist=155.5 / 0.705): + def __init__(self, nsim, paths, maxdist=155.5 / 0.705, + minmass=('M', 1e12)): self.nsim = nsim self.paths = paths # Read in the mmain catalogue of summed substructure @@ -417,4 +371,5 @@ class HaloCatalogue(BaseCatalogue): data = self.box.convert_from_boxunits(data, ['x', 'y', 'z', 'M']) mask = numpy.sqrt(data['x']**2 + data['y']**2 + data['z']**2) < maxdist + mask &= data[minmass[0]] > minmass[1] self._data = data[mask] diff --git a/csiborgtools/read/outsim.py b/csiborgtools/read/outsim.py deleted file mode 100644 index 91521fd..0000000 --- a/csiborgtools/read/outsim.py +++ /dev/null @@ -1,112 +0,0 @@ -# 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. -""" -I/O functions for analysing the CSiBORG realisations. -""" -from os import remove -from os.path import join - -import numpy -from tqdm import trange - - -def dump_split(arr, nsplit, nsnap, nsim, paths): - """ - Dump an array from a split. - - Parameters - ---------- - arr : n-dimensional or structured array - Array to be saved. - nsplit : int - Split index. - nsnap : int - Snapshot index. - nsim : int - IC realisation index. - paths : py:class`csiborgtools.read.CSiBORGPaths` - CSiBORG paths-handling object with set `n_sim` and `n_snap`. - - Returns - ------- - None - """ - fname = join(paths.temp_dumpdir, "ramses_out_{}_{}_{}.npy" - .format(str(nsim).zfill(5), str(nsnap).zfill(5), nsplit)) - numpy.save(fname, arr) - - -def combine_splits(nsplits, nsnap, nsim, part_reader, cols_add, - remove_splits=False, verbose=True): - """ - Combine results of many splits saved from `dump_split`. Identifies to which - clump the clumps in the split correspond to by matching their index. - Returns an array that contains the original clump data along with the newly - calculated quantities. - - Paramaters - ---------- - nsplits : int - Total number of clump splits. - nsnap : int - Snapshot index. - nsim : int - IC realisation index. - part_reader : py:class`csiborgtools.read.ParticleReadear` - CSiBORG particle reader. - cols_add : list of `(str, dtype)` - Colums to add. Must be formatted as, for example, - `[("npart", numpy.float64), ("totpartmass", numpy.float64)]`. - remove_splits : bool, optional - Whether to remove the splits files. By default `False`. - verbose : bool, optional - Verbosity flag. By default `True`. - - Returns - ------- - out : structured array - Clump array with appended results from the splits. - """ - clumps = part_reader.read_clumps(nsnap, nsim, cols=None) - # Get the old + new dtypes and create an empty array - descr = clumps.dtype.descr + cols_add - out = numpy.full(clumps.size, numpy.nan, dtype=descr) - for par in clumps.dtype.names: # Now put the old values into the array - out[par] = clumps[par] - - # Filename of splits data - froot = "ramses_out_{}_{}".format(str(nsim).zfill(5), str(nsnap).zfill(5)) - fname = join(part_reader.paths.temp_dumpdir, froot + "_{}.npy") - - # Iterate over splits and add to the output array - cols_add_names = [col[0] for col in cols_add] - iters = trange(nsplits) if verbose else range(nsplits) - for n in iters: - fnamesplit = fname.format(n) - arr = numpy.load(fnamesplit) - - # Check that all halo indices from the split are in the clump file - if not numpy.alltrue(numpy.isin(arr["index"], out["index"])): - raise KeyError("....") - # Mask of where to put the values from the split - mask = numpy.isin(out["index"], arr["index"]) - for par in cols_add_names: - out[par][mask] = arr[par] - - # Now remove this split - if remove_splits: - remove(fnamesplit) - - return out diff --git a/csiborgtools/read/readsim.py b/csiborgtools/read/readsim.py index 4035089..b709353 100644 --- a/csiborgtools/read/readsim.py +++ b/csiborgtools/read/readsim.py @@ -56,7 +56,7 @@ class ParticleReader: @paths.setter def paths(self, paths): - # assert isinstance(paths, CSiBORGPaths) # REMOVE + assert isinstance(paths, CSiBORGPaths) self._paths = paths def read_info(self, nsnap, nsim): @@ -87,7 +87,8 @@ class ParticleReader: keys = info[eqs - 1] vals = info[eqs + 1] - return {key: val for key, val in zip(keys, vals, strict=True)} + # trunk-ignore(ruff/B905) + return {key: val for key, val in zip(keys, vals)} def open_particle(self, nsnap, nsim, verbose=True): """ @@ -245,7 +246,8 @@ class ParticleReader: for cpu in iters: i = start_ind[cpu] j = nparts[cpu] - for (fname, fdtype) in zip(fnames, fdtypes, strict=True): + # trunk-ignore(ruff/B905) + for (fname, fdtype) in zip(fnames, fdtypes): if fname in pars_extract: out[fname][i:i + j] = self.read_sp(fdtype, partfiles[cpu]) else: diff --git a/csiborgtools/read/utils.py b/csiborgtools/read/utils.py index 1a7e140..d64db77 100644 --- a/csiborgtools/read/utils.py +++ b/csiborgtools/read/utils.py @@ -22,53 +22,59 @@ import numpy ############################################################################### -def cartesian_to_radec(x, y, z): +def cartesian_to_radec(X, indeg=True): """ - Calculate the radial distance, right ascension in [0, 360) degrees and - declination [-90, 90] degrees. Note, the observer should be placed in the - middle of the box. + Calculate the radial distance, RA, dec from Cartesian coordinates. Note, + RA is in range [0, 360) degrees and dec in range [-90, 90] degrees. Parameters ---------- - x, y, z : 1-dimensional arrays + X : 2-dimensional array `(nsamples, 3)` Cartesian coordinates. + indeg : bool, optional + Whether to return RA and DEC in degrees. Returns ------- - dist, ra, dec : 1-dimensional arrays - Radial distance, right ascension and declination. + out : 2-dimensional array `(nsamples, 3)` + Radial distance, RA and dec. """ - dist = numpy.sqrt(x**2 + y**2 + z**2) - dec = numpy.rad2deg(numpy.arcsin(z/dist)) - ra = numpy.rad2deg(numpy.arctan2(y, x)) - # Make sure RA in the correct range - ra[ra < 0] += 360 - return dist, ra, dec + x, y, z = X[:, 0], X[:, 1], X[:, 2] + dist = numpy.linalg.norm(X, axis=1) + dec = numpy.arcsin(z/dist) + ra = numpy.arctan2(y, x) + ra[ra < 0] += 2 * numpy.pi # Wrap RA to [0, 2pi) + if indeg: + ra = numpy.rad2deg(ra) + dec = numpy.rad2deg(dec) + return numpy.vstack([dist, ra, dec]).T -def radec_to_cartesian(dist, ra, dec, isdeg=True): +def radec_to_cartesian(X, isdeg=True): """ - Convert distance, right ascension and declination to Cartesian coordinates. + Calculate Cartesian coordinates from radial distance, RA, dec. Note, RA is + expected in range [0, 360) degrees and dec in range [-90, 90] degrees. Parameters ---------- - dist, ra, dec : 1-dimensional arrays - Spherical coordinates. + X : 2-dimensional array `(nsamples, 3)` + Radial distance, RA and dec. isdeg : bool, optional - Whether `ra` and `dec` are in degres. By default `True`. + Whether to return RA and DEC in degrees. Returns ------- - x, y, z : 1-dimensional arrays + out : 2-dimensional array `(nsamples, 3)` Cartesian coordinates. """ + dist, ra, dec = X[:, 0], X[:, 1], X[:, 2] if isdeg: ra = numpy.deg2rad(ra) dec = numpy.deg2rad(dec) - x = dist * numpy.cos(dec) * numpy.cos(ra) - y = dist * numpy.cos(dec) * numpy.sin(ra) - z = dist * numpy.sin(dec) - return x, y, z + x = numpy.cos(dec) * numpy.cos(ra) + y = numpy.cos(dec) * numpy.sin(ra) + z = numpy.sin(dec) + return dist * numpy.vstack([x, y, z]).T ############################################################################### diff --git a/notebooks/perseus.ipynb b/notebooks/perseus.ipynb new file mode 100644 index 0000000..9888e4e --- /dev/null +++ b/notebooks/perseus.ipynb @@ -0,0 +1,165 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "id": "5a38ed25", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-12T14:25:46.519408Z", + "start_time": "2023-04-12T14:25:03.003304Z" + }, + "scrolled": true + }, + "outputs": [], + "source": [ + "import sys\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import scienceplots\n", + "import astroquery\n", + "from tqdm import trange, tqdm\n", + "\n", + "sys.path.append(\"../\")\n", + "import csiborgtools\n", + "\n", + "%matplotlib widget \n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "79919270", + "metadata": {}, + "outputs": [], + "source": [ + "# # Norma\n", + "cluster = {\"RA\": (16 + 15 / 60 + 32.8 / 60**2) * 15,\n", + " \"DEC\": -60 + 54 / 60 + 30 / 60**2,\n", + " \"DIST\": 67.8}\n", + "\n", + "Xclust = np.array([cluster[\"DIST\"], cluster[\"RA\"], cluster[\"DEC\"]]).reshape(1, -1)" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "96740b90", + "metadata": {}, + "outputs": [], + "source": [ + "paths = csiborgtools.read.CSiBORGPaths(**csiborgtools.paths_glamdring)\n", + "nsims = paths.get_ics(False)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "90033fb1", + "metadata": {}, + "outputs": [], + "source": [ + "Xclust = np.array([cluster[\"DIST\"], cluster[\"RA\"], cluster[\"DEC\"]]).reshape(1, -1)" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "194baa83", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 101/101 [00:44<00:00, 2.25it/s]\n" + ] + } + ], + "source": [ + "matches = np.full(len(nsims), np.nan)\n", + "\n", + "for ii in trange(101):\n", + " cat = csiborgtools.read.HaloCatalogue(nsims[ii], paths, minmass=('M', 1e13))\n", + " dist, ind = cat.angular_neighbours(Xclust, ang_radius=5, rad_tolerance=10)\n", + " dist = dist[0]\n", + " ind = ind[0]\n", + "\n", + " if ind.size > 0:\n", + " matches[ii] = np.max(cat['M'][ind])\n" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "392f6eee", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "6ff7d5b9dc2f4b3fa6f563ee91f59655", + "version_major": 2, + "version_minor": 0 + }, + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAflklEQVR4nO3df3TV9X348Vf4kcCEBBMlISMBWltjtbAOK6a6ddXMyOH4o6Q/7HGWUk532hOpkLNV6dZST6fB7kxoK2D1ODw7p8zWP7CjntrDsonHs4AYxo6uk1ZrD1hMqG4kiCMw8vn+0a+p0UT5kZub3Pfjcc49x3zuzcfXy+Anz3NzbyjKsiwLAACSMS7fAwAAMLIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYgQgAEBiBCAAQGIEIABAYibke4CxrK+vLw4cOBBTp06NoqKifI8DAJyELMvi8OHDUV1dHePGpflcmAA8AwcOHIiampp8jwEAnIb9+/fHzJkz8z1GXgjAMzB16tSI+O0foNLS0jxPAwCcjJ6enqipqen/Pp4iAXgG3vixb2lpqQAEgDEm5ZdvpfmDbwCAhAlAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDETMj3AABA7s2+7dF8j3DKfrVmUb5HKFieAQQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEhMwQbgN77xjSgqKhpwq6ur67//6NGj0dzcHBUVFTFlypRoamqKrq6uPE4MADAyCjYAIyIuvPDCePnll/tvTz75ZP99K1eujK1bt8bDDz8c27dvjwMHDsTixYvzOC0AwMiYkO8BcmnChAlRVVX1tuPd3d3xwAMPxObNm+OKK66IiIhNmzbFBRdcEDt27IhLL710pEcFABgxBf0M4C9+8Yuorq6O97znPXHjjTfGvn37IiKio6Mjjh8/Hg0NDf2Prauri9ra2mhvbx/yfL29vdHT0zPgBgAw1hRsAC5YsCAefPDBeOyxx2Ljxo3x4osvxh/90R/F4cOHo7OzM4qLi2PatGkDPqeysjI6OzuHPGdra2uUlZX132pqanK8BQDA8CvYHwEvXLiw/5/nzp0bCxYsiFmzZsUPf/jDmDx58mmdc9WqVdHS0tL/cU9PjwgEAMacgn0G8K2mTZsW73//++P555+PqqqqOHbsWBw6dGjAY7q6ugZ9zeAbSkpKorS0dMANAGCsSSYAX3vttXjhhRdixowZMX/+/Jg4cWK0tbX13793797Yt29f1NfX53FKAIDcK9gfAf/FX/xFXHPNNTFr1qw4cOBArF69OsaPHx+f+cxnoqysLJYtWxYtLS1RXl4epaWlsXz58qivr/cOYACg4BVsAL700kvxmc98Jl599dU499xz4/LLL48dO3bEueeeGxERa9eujXHjxkVTU1P09vZGY2NjbNiwIc9TAwDkXlGWZVm+hxirenp6oqysLLq7u70eEIBRbfZtj+Z7hFP2qzWLcnJe378Teg0gAAC/JQABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABIjAAEAEiMAAQASIwABABKTRACuWbMmioqKYsWKFf3Hjh49Gs3NzVFRURFTpkyJpqam6Orqyt+QAAAjpOADcNeuXfG9730v5s6dO+D4ypUrY+vWrfHwww/H9u3b48CBA7F48eI8TQkAMHIKOgBfe+21uPHGG+P++++Ps88+u/94d3d3PPDAA3H33XfHFVdcEfPnz49NmzbFv/3bv8WOHTvyODEAQO4VdAA2NzfHokWLoqGhYcDxjo6OOH78+IDjdXV1UVtbG+3t7UOer7e3N3p6egbcAADGmgn5HiBXHnroodi9e3fs2rXrbfd1dnZGcXFxTJs2bcDxysrK6OzsHPKcra2tcfvttw/3qAAAI6ognwHcv39/3HLLLfH9738/Jk2aNGznXbVqVXR3d/ff9u/fP2znBgAYKQUZgB0dHXHw4MH4wz/8w5gwYUJMmDAhtm/fHt/5zndiwoQJUVlZGceOHYtDhw4N+Lyurq6oqqoa8rwlJSVRWlo64AYAMNYU5I+Ar7zyynjmmWcGHFu6dGnU1dXFrbfeGjU1NTFx4sRoa2uLpqamiIjYu3dv7Nu3L+rr6/MxMgDAiCnIAJw6dWpcdNFFA46dddZZUVFR0X982bJl0dLSEuXl5VFaWhrLly+P+vr6uPTSS/MxMgDAiCnIADwZa9eujXHjxkVTU1P09vZGY2NjbNiwId9jAQDkXFGWZVm+hxirenp6oqysLLq7u70eEIBRbfZtj+Z7hFP2qzWLcnJe378L9E0gAAAMTQACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACRGAAIAJEYAAgAkRgACACSmYANw48aNMXfu3CgtLY3S0tKor6+Pn/zkJ/33Hz16NJqbm6OioiKmTJkSTU1N0dXVlceJAQBGRsEG4MyZM2PNmjXR0dERTz/9dFxxxRVx3XXXxX/+539GRMTKlStj69at8fDDD8f27dvjwIEDsXjx4jxPDQCQe0VZlmX5HmKklJeXx9/+7d/GJz7xiTj33HNj8+bN8YlPfCIiIp577rm44IILor29PS699NKTOl9PT0+UlZVFd3d3lJaW5nJ0ADgjs297NN8jnLJfrVmUk/P6/l3AzwC+2YkTJ+Khhx6KI0eORH19fXR0dMTx48ejoaGh/zF1dXVRW1sb7e3tQ56nt7c3enp6BtwAAMaagg7AZ555JqZMmRIlJSXxxS9+MbZs2RIf+MAHorOzM4qLi2PatGkDHl9ZWRmdnZ1Dnq+1tTXKysr6bzU1NTneAABg+BV0AJ5//vmxZ8+e2LlzZ3zpS1+KJUuWxM9+9rPTPt+qVauiu7u7/7Z///5hnBYAYGRMyPcAuVRcXBznnXdeRETMnz8/du3aFd/+9rfj05/+dBw7diwOHTo04FnArq6uqKqqGvJ8JSUlUVJSkuuxAQByqqCfAXyrvr6+6O3tjfnz58fEiROjra2t/769e/fGvn37or6+Po8TAgDkXsE+A7hq1apYuHBh1NbWxuHDh2Pz5s3x+OOPx09/+tMoKyuLZcuWRUtLS5SXl0dpaWksX7486uvrT/odwAAAY1XBBuDBgwfjs5/9bLz88stRVlYWc+fOjZ/+9Kfxp3/6pxERsXbt2hg3blw0NTVFb29vNDY2xoYNG/I8NQBA7iX1ewCHm98jBMBY4fcA/o7v34m9BhAAAAEIAJAcAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJCYCfkeAADGmtm3PZrvEeCMeAYQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDEFG4Ctra3x4Q9/OKZOnRrTp0+P66+/Pvbu3TvgMUePHo3m5uaoqKiIKVOmRFNTU3R1deVpYgCAkVGwAbh9+/Zobm6OHTt2xLZt2+L48eNx1VVXxZEjR/ofs3Llyti6dWs8/PDDsX379jhw4EAsXrw4j1MDAOTehHwPkCuPPfbYgI8ffPDBmD59enR0dMQf//EfR3d3dzzwwAOxefPmuOKKKyIiYtOmTXHBBRfEjh074tJLL83H2AAAOVewzwC+VXd3d0RElJeXR0RER0dHHD9+PBoaGvofU1dXF7W1tdHe3j7oOXp7e6Onp2fADQBgrEkiAPv6+mLFihVx2WWXxUUXXRQREZ2dnVFcXBzTpk0b8NjKysro7Owc9Dytra1RVlbWf6upqcn16AAAwy6JAGxubo5nn302HnrooTM6z6pVq6K7u7v/tn///mGaEABg5BTsawDfcPPNN8ePf/zjeOKJJ2LmzJn9x6uqquLYsWNx6NChAc8CdnV1RVVV1aDnKikpiZKSklyPDACQUwX7DGCWZXHzzTfHli1b4l/+5V9izpw5A+6fP39+TJw4Mdra2vqP7d27N/bt2xf19fUjPS4AwIgp2GcAm5ubY/PmzfGjH/0opk6d2v+6vrKyspg8eXKUlZXFsmXLoqWlJcrLy6O0tDSWL18e9fX13gEMABS0gg3AjRs3RkTEn/zJnww4vmnTpvjc5z4XERFr166NcePGRVNTU/T29kZjY2Ns2LBhhCcFABhZBRuAWZa962MmTZoU69evj/Xr14/ARAAAo0PBvgYQAIDBCUAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQUbAA+8cQTcc0110R1dXUUFRXFI488MuD+LMvi61//esyYMSMmT54cDQ0N8Ytf/CI/wwIAjKCCDcAjR47EvHnzYv369YPe/61vfSu+853vxL333hs7d+6Ms846KxobG+Po0aMjPCkAwMiakO8BcmXhwoWxcOHCQe/LsizWrVsXf/3Xfx3XXXddRET8wz/8Q1RWVsYjjzwSN9xww0iOCgAwogr2GcB38uKLL0ZnZ2c0NDT0HysrK4sFCxZEe3v7kJ/X29sbPT09A24AAGNNkgHY2dkZERGVlZUDjldWVvbfN5jW1tYoKyvrv9XU1OR0TgCAXEgyAE/XqlWroru7u/+2f//+fI8EAHDKkgzAqqqqiIjo6uoacLyrq6v/vsGUlJREaWnpgBsAwFiTZADOmTMnqqqqoq2trf9YT09P7Ny5M+rr6/M4GQBA7hXsu4Bfe+21eP755/s/fvHFF2PPnj1RXl4etbW1sWLFivibv/mbeN/73hdz5syJr33ta1FdXR3XX399/oYGABgBBRuATz/9dHzsYx/r/7ilpSUiIpYsWRIPPvhgfOUrX4kjR47En//5n8ehQ4fi8ssvj8ceeywmTZqUr5EBAEZEUZZlWb6HGKt6enqirKwsuru7vR4QICGzb3s03yMk4VdrFuXkvL5/J/oaQACAlAlAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxEzI9wAApG32bY/mewRIjmcAAQASIwABABIjAAEAEiMAAQASIwABABLjXcAAQ/DuVKBQeQYQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMQIQACAxAhAAIDECEAAgMRPyPQCQhtm3PZrvEQD4/zwDCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkBgBCACQGAEIAJAYAQgAkJgJ+R6Aoc2+7dF8j5CEX61ZlO8RTpk/GwCcCc8AAgAkJvkAXL9+fcyePTsmTZoUCxYsiKeeeirfIwEA5FTSAfiDH/wgWlpaYvXq1bF79+6YN29eNDY2xsGDB/M9GgBAziQdgHfffXd84QtfiKVLl8YHPvCBuPfee+P3fu/34u///u/zPRoAQM4k+yaQY8eORUdHR6xatar/2Lhx46KhoSHa29sH/Zze3t7o7e3t/7i7uzsiInp6enIyY1/v6zk5LwPl6uuXS/5sACnI1fX5jfNmWZaT848FyQbgK6+8EidOnIjKysoBxysrK+O5554b9HNaW1vj9ttvf9vxmpqanMzIyChbl+8JABhMrq/Phw8fjrKystz+S0apZAPwdKxatSpaWlr6P+7r64v//u//joqKiigqKsrjZMOvp6cnampqYv/+/VFaWprvcYZdIe9XyLtFFPZ+hbxbRGHvV8i7RRTeflmWxeHDh6O6ujrfo+RNsgF4zjnnxPjx46Orq2vA8a6urqiqqhr0c0pKSqKkpGTAsWnTpuVqxFGhtLS0IP5nH0oh71fIu0UU9n6FvFtEYe9XyLtFFNZ+qT7z94Zk3wRSXFwc8+fPj7a2tv5jfX190dbWFvX19XmcDAAgt5J9BjAioqWlJZYsWRIXX3xxXHLJJbFu3bo4cuRILF26NN+jAQDkTNIB+OlPfzp+85vfxNe//vXo7OyMP/iDP4jHHnvsbW8MSVFJSUmsXr36bT/yLhSFvF8h7xZR2PsV8m4Rhb1fIe8WUfj7pagoS/k90AAACUr2NYAAAKkSgAAAiRGAAACJEYAAAIkRgAl44okn4pprronq6uooKiqKRx55ZMD93/jGN6Kuri7OOuusOPvss6OhoSF27tz5juc8ceJEfO1rX4s5c+bE5MmT473vfW9885vfzMvfq/hu+73ZF7/4xSgqKop169a963nXr18fs2fPjkmTJsWCBQviqaeeGr6hT1IudmttbY0Pf/jDMXXq1Jg+fXpcf/31sXfv3uEd/CTl6mv3hjVr1kRRUVGsWLHijGc9Vbna7de//nX82Z/9WVRUVMTkyZPjgx/8YDz99NPDN/hJysV+o+W68m67fe5zn4uioqIBt6uvvvpdzzsarikRudlvNF1XODkCMAFHjhyJefPmxfr16we9//3vf3/cc8898cwzz8STTz4Zs2fPjquuuip+85vfDHnOu+66KzZu3Bj33HNP/Nd//Vfcdddd8a1vfSu++93v5mqNIb3bfm/YsmVL7Nix46T+6p8f/OAH0dLSEqtXr47du3fHvHnzorGxMQ4ePDhcY5+UXOy2ffv2aG5ujh07dsS2bdvi+PHjcdVVV8WRI0eGa+yTlov93rBr16743ve+F3Pnzj3TMU9LLnb7n//5n7jsssti4sSJ8ZOf/CR+9rOfxd/93d/F2WefPVxjn7Rc7Ddarisns9vVV18dL7/8cv/tH//xH9/xnKPlmhKRm/1G03WFk5SRlIjItmzZ8o6P6e7uziIi++d//uchH7No0aLs85///IBjixcvzm688cbhGPO0DbXfSy+9lP3+7/9+9uyzz2azZs3K1q5d+47nueSSS7Lm5ub+j0+cOJFVV1dnra2twzzxyRuu3d7q4MGDWURk27dvH55BT9Nw7nf48OHsfe97X7Zt27bsox/9aHbLLbcM+7ynYrh2u/XWW7PLL788N0OegeHabzReVwbbbcmSJdl11113SucZjdeULBu+/d5qtFxXGJpnABng2LFjcd9990VZWVnMmzdvyMd95CMfiba2tvj5z38eERH/8R//EU8++WQsXLhwpEY9aX19fXHTTTfFX/7lX8aFF174ro8/duxYdHR0RENDQ/+xcePGRUNDQ7S3t+dy1FN2qrsNpru7OyIiysvLh3O0YXG6+zU3N8eiRYsGfA1Hm9PZ7Z/+6Z/i4osvjk9+8pMxffr0+NCHPhT3339/jic9Paez31i6rjz++OMxffr0OP/88+NLX/pSvPrqq0M+dixdU95wKvsNZjRfV/itpP8mEH7nxz/+cdxwww3x+uuvx4wZM2Lbtm1xzjnnDPn42267LXp6eqKuri7Gjx8fJ06ciDvuuCNuvPHGEZz65Nx1110xYcKE+PKXv3xSj3/llVfixIkTb/sbYSorK+O5557LxYin7VR3e6u+vr5YsWJFXHbZZXHRRRcN83Rn7nT2e+ihh2L37t2xa9euHE525k5nt1/+8pexcePGaGlpia9+9auxa9eu+PKXvxzFxcWxZMmSHE576k5nv7FyXbn66qtj8eLFMWfOnHjhhRfiq1/9aixcuDDa29tj/Pjxb3v8WLqmRJz6fm812q8r/JYAJCIiPvaxj8WePXvilVdeifvvvz8+9alPxc6dO2P69OmDPv6HP/xhfP/734/NmzfHhRdeGHv27IkVK1ZEdXX1qPpG1NHREd/+9rdj9+7dUVRUlO9xhtVw7Nbc3BzPPvtsPPnkk8M83Zk7nf32798ft9xyS2zbti0mTZqU4wlP3+l+7fr6+uLiiy+OO++8MyIiPvShD8Wzzz4b9957b0H8fzdWris33HBD/z9/8IMfjLlz58Z73/veePzxx+PKK6/M42TD40z3G83XFd4k3z+DZmTFSbwGMMuy7LzzzsvuvPPOIe+fOXNmds899ww49s1vfjM7//zzz3TEM/LW/dauXZsVFRVl48eP779FRDZu3Lhs1qxZg56jt7c3Gz9+/Nv+O332s5/Nrr322twN/y6GY7c3a25uzmbOnJn98pe/zN3Qp2A49tuyZUsWEW/7nDfO83//938js8xbDNfXrra2Nlu2bNmAYxs2bMiqq6tzNPnJGa79RuN15WSvmeecc0527733DnrfaL2mZNnw7Pdmo+26wtA8A8ig+vr6ore3d8j7X3/99Rg3buBLSMePHx99fX25Hu2U3HTTTW97HVhjY2PcdNNNsXTp0kE/p7i4OObPnx9tbW1x/fXXR8Rv/3u0tbXFzTffnOuRT9rp7BYRkWVZLF++PLZs2RKPP/54zJkzJ9ejnpbT2e/KK6+MZ555ZsCxpUuXRl1dXdx6660n9eOrkXC6X7vLLrvsbb9a4+c//3nMmjUrJ3OertPdb6xcV97qpZdeildffTVmzJgx6P1j5ZoylHfbL2LsXFf4HQGYgNdeey2ef/75/o9ffPHF2LNnT5SXl0dFRUXccccdce2118aMGTPilVdeifXr18evf/3r+OQnP9n/OVdeeWV8/OMf779YXXPNNXHHHXdEbW1tXHjhhfHv//7vcffdd8fnP//5UbVfbW1tVFRUDHj8xIkTo6qqKs4///z+Y2/dr6WlJZYsWRIXX3xxXHLJJbFu3bo4cuTIO37zyoVc7Nbc3BybN2+OH/3oRzF16tTo7OyMiIiysrKYPHnyCGz1O8O939SpU9/2mqOzzjorKioqRvy1SLn42q1cuTI+8pGPxJ133hmf+tSn4qmnnor77rsv7rvvvpFZ6k1ysd9oua68027l5eVx++23R1NTU1RVVcULL7wQX/nKV+K8886LxsbGIXcbLdeUXO03mq4rnKR8PwVJ7v3rv/5rFhFvuy1ZsiT73//93+zjH/94Vl1dnRUXF2czZszIrr322uypp54acI5Zs2Zlq1ev7v+4p6cnu+WWW7La2tps0qRJ2Xve857sr/7qr7Le3t4R3u6d9xvMYL+O4q37ZVmWffe7381qa2uz4uLi7JJLLsl27NiRmwXeQS52G+x8EZFt2rQpZ3sMJVdfuzfL16+BydVuW7duzS666KKspKQkq6ury+67777cLPAucrHfaLmuvNNur7/+enbVVVdl5557bjZx4sRs1qxZ2Re+8IWss7PzHXfLstFxTcmy3Ow3mq4rnJyiLMvDX90AAEDe+D2AAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJEYAAAIkRgAAAiRGAAACJ+X/SiIf6vFeUOQAAAABJRU5ErkJggg==", + "text/html": [ + "\n", + "
\n", + "
\n", + " Figure\n", + "
\n", + " \n", + "
\n", + " " + ], + "text/plain": [ + "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x = np.log10(matches[~np.isnan(matches)])\n", + "\n", + "\n", + "plt.figure()\n", + "plt.hist(x, bins=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9cc660b0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "venv_galomatch", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}