diff --git a/.gitignore b/.gitignore index 1814f26..de42f81 100644 --- a/.gitignore +++ b/.gitignore @@ -15,4 +15,4 @@ build/* csiborgtools.egg-info/* Pylians3/* scripts/plot_correlation.ipynb -scripts/python.sh +scripts/*.sh diff --git a/README.md b/README.md index 469d4a9..88587e5 100644 --- a/README.md +++ b/README.md @@ -7,12 +7,20 @@ ## Project Clustering -- [ ] Add uncertainty to the kNN-CDF autocorrelation? -- [ ] Add kNN-CDF differences. -- [ ] Add reading halo catalogues at higher redshifts. -- [x] Add the joint kNN-CDF calculation. -- [x] Make kNN-CDF more memory friendly if generating many randoms. +### Longterm +- [ ] Add uncertainty to the kNN-CDF autocorrelation? +- [ ] Add reading halo catalogues at higher redshifts. + + +### April 9 2023 Sunday +- [x] Add normalised marks calculation. +- [x] Add normalised marks to the submission scripts. +- [x] Verify analytical formula for the kNN of a uniform field. +- [x] For the cross-correlation try making the second field randoms. +- [ ] Clean up the reader code. +- [x] Correct the crossing script. +- [ ] Get started with the 2PCF calculation. ## Project Environmental Dependence - [ ] Add gradient and Hessian of the overdensity field. diff --git a/csiborgtools/__init__.py b/csiborgtools/__init__.py index 5c2486d..1ff0324 100644 --- a/csiborgtools/__init__.py +++ b/csiborgtools/__init__.py @@ -12,4 +12,4 @@ # 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. -from csiborgtools import (read, match, utils, units, fits, field) # noqa +from csiborgtools import (read, match, utils, units, fits, field, clustering) # noqa diff --git a/csiborgtools/match/correlation.py b/csiborgtools/clustering/2pcf.py similarity index 69% rename from csiborgtools/match/correlation.py rename to csiborgtools/clustering/2pcf.py index 5c376a9..16e0284 100644 --- a/csiborgtools/match/correlation.py +++ b/csiborgtools/clustering/2pcf.py @@ -1,4 +1,4 @@ -# Copyright (C) 2022 Richard Stiskalek +# Copyright (C) 2023 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 @@ -12,58 +12,15 @@ # 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. +""" +2PCF calculation. + +NOTE: This is an old script that needs to be updated. +""" import numpy from Corrfunc.mocks import DDtheta_mocks from Corrfunc.utils import convert_3d_counts_to_cf -from warnings import warn - - -def get_randoms_sphere(N, seed=42): - """ - Generate random points on a sphere. - - Parameters - ---------- - N : int - Number of points. - seed : int - Random seed. - - Returns - ------- - ra : 1-dimensional array - Right ascension in :math:`[0, 360)` degrees. - dec : 1-dimensional array - Declination in :math:`[-90, 90]` degrees. - """ - gen = numpy.random.default_rng(seed) - ra = gen.random(N) * 360 - dec = numpy.rad2deg(numpy.arcsin(2 * (gen.random(N) - 0.5))) - return ra, dec - - -def wrapRA(ra, degrees=True): - """ - Wrap the right ascension from :math:`[-180, 180)` to :math`[0, 360)` - degrees or equivalently if `degrees=False` in radians. - - Paramaters - ---------- - ra : 1-dimensional array - Right ascension values. - degrees : float, optional - Whether the right ascension is in degrees. - - Returns - ------- - ra : 1-dimensional array - Wrapped around right ascension. - """ - mask = ra < 0 - if numpy.sum(mask) == 0: - warn("No negative right ascension found.") - ra[mask] += 360 if degrees else 2 * numpy.pi - return ra +from .utils import (rvs_on_sphere, wrapRA) def sphere_angular_tpcf(bins, RA1, DEC1, RA2=None, DEC2=None, nthreads=1, @@ -113,11 +70,11 @@ def sphere_angular_tpcf(bins, RA1, DEC1, RA2=None, DEC2=None, nthreads=1, NR1 = ND1 * Nmult NR2 = ND2 * Nmult # Generate randoms. Note that these are over the sphere! - randRA1, randDEC1 = get_randoms_sphere(NR1, seed1) - randRA2, randDEC2 = get_randoms_sphere(NR2, seed2) + randRA1, randDEC1 = rvs_on_sphere(NR1, indeg=True, random_state=seed1) + randRA2, randDEC2 = rvs_on_sphere(NR2, indeg=True, random_state=seed2) # Wrap RA - RA1 = wrapRA(numpy.copy(RA1)) - RA2 = wrapRA(numpy.copy(RA2)) + RA1 = wrapRA(numpy.copy(RA1), indeg=True) + RA2 = wrapRA(numpy.copy(RA2), indeg=True) # Calculate pairs D1D2 = DDtheta_mocks(0, nthreads, bins, RA1, DEC1, RA2=RA2, DEC2=DEC2) D1R2 = DDtheta_mocks(0, nthreads, bins, RA1, DEC1, @@ -127,4 +84,4 @@ def sphere_angular_tpcf(bins, RA1, DEC1, RA2=None, DEC2=None, nthreads=1, R1R2 = DDtheta_mocks(0, nthreads, bins, randRA1, randDEC1, RA2=randRA2, DEC2=randDEC2) # Convert to the CF - return convert_3d_counts_to_cf(ND1, ND2, NR1, NR2, D1D2, D1R2, D2R1, R1R2) + return convert_3d_counts_to_cf(ND1, ND2, NR1, NR2, D1D2, D1R2, D2R1, R1R2) \ No newline at end of file diff --git a/csiborgtools/clustering/__init__.py b/csiborgtools/clustering/__init__.py new file mode 100644 index 0000000..49b32f6 --- /dev/null +++ b/csiborgtools/clustering/__init__.py @@ -0,0 +1,16 @@ +# 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. +from .knn import kNN_CDF # noqa +from .utils import (RVSinsphere, RVSinbox, RVSonsphere, BaseRVS, normalised_marks) # noqa diff --git a/csiborgtools/match/knn.py b/csiborgtools/clustering/knn.py similarity index 70% rename from csiborgtools/match/knn.py rename to csiborgtools/clustering/knn.py index 8904d44..e57beb7 100644 --- a/csiborgtools/match/knn.py +++ b/csiborgtools/clustering/knn.py @@ -18,52 +18,16 @@ kNN-CDF calculation import numpy from scipy.interpolate import interp1d from scipy.stats import binned_statistic -from tqdm import tqdm +from .utils import BaseRVS class kNN_CDF: - """ - Object to calculate the kNN-CDF for a set of CSiBORG halo catalogues from - their kNN objects. - """ - @staticmethod - def rvs_in_sphere(nsamples, R, random_state=42, dtype=numpy.float32): - """ - Generate random samples in a sphere of radius `R` centered at the - origin. - - Parameters - ---------- - nsamples : int - Number of samples to generate. - R : float - Radius of the sphere. - random_state : int, optional - Random state for the random number generator. - dtype : numpy dtype, optional - Data type, by default `numpy.float32`. - - Returns - ------- - samples : 2-dimensional array of shape `(nsamples, 3)` - """ - gen = numpy.random.default_rng(random_state) - # Sample spherical coordinates - r = gen.uniform(0, 1, nsamples).astype(dtype)**(1/3) * R - theta = 2 * numpy.arcsin(gen.uniform(0, 1, nsamples).astype(dtype)) - phi = 2 * numpy.pi * gen.uniform(0, 1, nsamples).astype(dtype) - # Convert to cartesian coordinates - x = r * numpy.sin(theta) * numpy.cos(phi) - y = r * numpy.sin(theta) * numpy.sin(phi) - z = r * numpy.cos(theta) - - return numpy.vstack([x, y, z]).T - + """Object to calculate the kNN-CDF statistic.""" @staticmethod def cdf_from_samples(r, rmin=None, rmax=None, neval=None, dtype=numpy.float32): """ - Calculate the CDF from samples. + Calculate the kNN-CDF from a sampled PDF. Parameters ---------- @@ -128,22 +92,21 @@ class kNN_CDF: corr[k, :] = joint_cdf[k, :] - cdf0[k, :] * cdf1[k, :] return corr - def brute_cdf(self, knn, nneighbours, Rmax, nsamples, rmin, rmax, neval, + def brute_cdf(self, knn, rvs_gen, nneighbours, nsamples, rmin, rmax, neval, random_state=42, dtype=numpy.float32): """ - Calculate the CDF for a kNN of CSiBORG halo catalogues without batch - sizing. This can become memory intense for large numbers of randoms - and, therefore, is only for testing purposes. + Calculate the kNN-CDF without batch sizing. This can become memory + intense for large numbers of randoms and, therefore, is primarily for + testing purposes. Parameters ---------- - knns : `sklearn.neighbors.NearestNeighbors` - kNN of CSiBORG halo catalogues. + knn : `sklearn.neighbors.NearestNeighbors` + Catalogue NN object. + rvs_gen : :py:class:`csiborgtools.clustering.BaseRVS` + Uniform RVS generator matching `knn`. neighbours : int Maximum number of neighbours to use for the kNN-CDF calculation. - Rmax : float - Maximum radius of the sphere in which to sample random points for - the knn-CDF calculation. This should match the CSiBORG catalogues. nsamples : int Number of random points to sample for the knn-CDF calculation. rmin : float @@ -164,7 +127,8 @@ class kNN_CDF: cdfs : 2-dimensional array CDFs evaluated at `rs`. """ - rand = self.rvs_in_sphere(nsamples, Rmax, random_state=random_state) + assert isinstance(rvs_gen, BaseRVS) + rand = rvs_gen(nsamples, random_state=random_state) dist, __ = knn.kneighbors(rand, nneighbours) dist = dist.astype(dtype) @@ -177,18 +141,20 @@ class kNN_CDF: cdf = numpy.asanyarray(cdf) return rs, cdf - def joint(self, knn0, knn1, nneighbours, Rmax, nsamples, rmin, rmax, + def joint(self, knn0, knn1, rvs_gen, nneighbours, nsamples, rmin, rmax, neval, batch_size=None, random_state=42, dtype=numpy.float32): """ - Calculate the joint CDF for two kNNs of CSiBORG halo catalogues. + Calculate the joint knn-CDF. Parameters ---------- knn0 : `sklearn.neighbors.NearestNeighbors` instance - kNN of the first CSiBORG halo catalogue. + NN object of the first catalogue. knn1 : `sklearn.neighbors.NearestNeighbors` instance - kNN of the second CSiBORG halo catalogue. + NN object of the second catalogue. + rvs_gen : :py:class:`csiborgtools.clustering.BaseRVS` + Uniform RVS generator matching `knn1` and `knn2`. neighbours : int Maximum number of neighbours to use for the kNN-CDF calculation. Rmax : float @@ -222,6 +188,7 @@ class kNN_CDF: joint_cdf : 2-dimensional array Joint CDF evaluated at `rs`. """ + assert isinstance(rvs_gen, BaseRVS) batch_size = nsamples if batch_size is None else batch_size assert nsamples >= batch_size nbatches = nsamples // batch_size @@ -233,8 +200,7 @@ class kNN_CDF: jointdist = numpy.zeros((batch_size, 2), dtype=dtype) for j in range(nbatches): - rand = self.rvs_in_sphere(batch_size, Rmax, - random_state=random_state + j) + rand = rvs_gen(batch_size, random_state=random_state + j) dist0, __ = knn0.kneighbors(rand, nneighbours) dist1, __ = knn1.kneighbors(rand, nneighbours) @@ -269,21 +235,19 @@ class kNN_CDF: rs = (bins[1:] + bins[:-1]) / 2 # Bin centers return rs, cdf0, cdf1, joint_cdf - def __call__(self, *knns, nneighbours, Rmax, nsamples, rmin, rmax, neval, - batch_size=None, verbose=True, random_state=42, - dtype=numpy.float32): + def __call__(self, knn, rvs_gen, nneighbours, nsamples, rmin, rmax, neval, + batch_size=None, random_state=42, dtype=numpy.float32): """ Calculate the CDF for a set of kNNs of CSiBORG halo catalogues. Parameters ---------- - *knns : `sklearn.neighbors.NearestNeighbors` instances - kNNs of CSiBORG halo catalogues. + knn : `sklearn.neighbors.NearestNeighbors` + Catalogue NN object. + rvs_gen : :py:class:`csiborgtools.clustering.BaseRVS` + Uniform RVS generator matching `knn1` and `knn2`. neighbours : int Maximum number of neighbours to use for the kNN-CDF calculation. - Rmax : float - Maximum radius of the sphere in which to sample random points for - the knn-CDF calculation. This should match the CSiBORG catalogues. nsamples : int Number of random points to sample for the knn-CDF calculation. rmin : float @@ -296,8 +260,6 @@ class kNN_CDF: Number of random points to sample in each batch. By default equal to `nsamples`, however recommeded to be smaller to avoid requesting too much memory, - verbose : bool, optional - Verbosity flag. random_state : int, optional Random state for the random number generator. dtype : numpy dtype, optional @@ -307,33 +269,30 @@ class kNN_CDF: ------- rs : 1-dimensional array Distances at which the CDF is evaluated. - cdfs : 2 or 3-dimensional array - CDFs evaluated at `rs`. + cdf : 2-dimensional array + CDF evaluated at `rs`. """ + assert isinstance(rvs_gen, BaseRVS) batch_size = nsamples if batch_size is None else batch_size assert nsamples >= batch_size nbatches = nsamples // batch_size # Preallocate the bins and the CDF array bins = numpy.logspace(numpy.log10(rmin), numpy.log10(rmax), neval) - cdfs = numpy.zeros((len(knns), nneighbours, neval - 1), dtype=dtype) - for i, knn in enumerate(tqdm(knns) if verbose else knns): - for j in range(nbatches): - rand = self.rvs_in_sphere(batch_size, Rmax, - random_state=random_state + j) - dist, __ = knn.kneighbors(rand, nneighbours) + cdf = numpy.zeros((nneighbours, neval - 1), dtype=dtype) + for i in range(nbatches): + rand = rvs_gen(batch_size, random_state=random_state + i) + dist, __ = knn.kneighbors(rand, nneighbours) - for k in range(nneighbours): # Count for each neighbour - _counts, __, __ = binned_statistic( - dist[:, k], dist[:, k], bins=bins, statistic="count", - range=(rmin, rmax)) - cdfs[i, k, :] += _counts + for k in range(nneighbours): # Count for each neighbour + _counts, __, __ = binned_statistic( + dist[:, k], dist[:, k], bins=bins, statistic="count", + range=(rmin, rmax)) + cdf[k, :] += _counts - cdfs = numpy.cumsum(cdfs, axis=-1) # Cumulative sum, i.e. the CDF - for i in range(len(knns)): - for k in range(nneighbours): - cdfs[i, k, :] /= cdfs[i, k, -1] + cdf = numpy.cumsum(cdf, axis=-1) # Cumulative sum, i.e. the CDF + for k in range(nneighbours): + cdf[k, :] /= cdf[k, -1] rs = (bins[1:] + bins[:-1]) / 2 # Bin centers - cdfs = cdfs[0, ...] if len(knns) == 1 else cdfs - return rs, cdfs + return rs, cdf diff --git a/csiborgtools/clustering/utils.py b/csiborgtools/clustering/utils.py new file mode 100644 index 0000000..839e6cc --- /dev/null +++ b/csiborgtools/clustering/utils.py @@ -0,0 +1,193 @@ +# 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. +"""Clustering support functions.""" +from abc import (ABC, abstractmethod) +from warnings import warn +import numpy + + +############################################################################### +# Random points # +############################################################################### + + +class BaseRVS(ABC): + """ + Base RVS generator. + """ + @abstractmethod + def __call__(self, nsamples, random_state, dtype): + """ + Generate RVS. + + Parameters + ---------- + nsamples : int + Number of samples to generate. + random_state : int, optional + Random state for the random number generator. + dtype : numpy dtype, optional + Data type, by default `numpy.float32`. + + Returns + ------- + samples : 2-dimensional array of shape `(nsamples, ndim)` + """ + pass + + +class RVSinsphere(BaseRVS): + """ + Generator of uniform RVS in a sphere of radius `R` in Cartesian + coordinates centered at the origin. + + Parameters + ---------- + R : float + Radius of the sphere. + """ + def __init__(self, R): + assert R > 0, "Radius must be positive." + self.R = R + BaseRVS.__init__(self) + + def __call__(self, nsamples, random_state=42, dtype=numpy.float32): + gen = numpy.random.default_rng(random_state) + # Spherical + r = gen.random(nsamples, dtype=dtype)**(1/3) * self.R + theta = 2 * numpy.arcsin(gen.random(nsamples, dtype=dtype)) + phi = 2 * numpy.pi * gen.random(nsamples, dtype=dtype) + # Cartesian + x = r * numpy.sin(theta) * numpy.cos(phi) + y = r * numpy.sin(theta) * numpy.sin(phi) + z = r * numpy.cos(theta) + return numpy.vstack([x, y, z]).T + + +class RVSinbox(BaseRVS): + """ + Generator of uniform RVS in a box of width `L` in Cartesian coordinates in + :math:`[0, L]^3`. + + Parameters + ---------- + width : float + Width of the box. + """ + def __init__(self, width): + assert width > 0, "Width must be positive." + self.width = width + BaseRVS.__init__(self) + + def __call__(self, nsamples, random_state=42, dtype=numpy.float32): + gen = numpy.random.default_rng(random_state) + x = gen.random(nsamples, dtype=dtype) + y = gen.random(nsamples, dtype=dtype) + z = gen.random(nsamples, dtype=dtype) + return self.width * numpy.vstack([x, y, z]).T + + +class RVSonsphere(BaseRVS): + """ + Generator of uniform RVS on the surface of a unit sphere. RA is in + :math:`[0, 2\pi)` and dec in :math:`[-\pi / 2, \pi / 2]`, respectively. + If `indeg` is `True` then converted to degrees. + + Parameters + ---------- + indeg : bool + Whether to generate the right ascension and declination in degrees. + """ + def __init__(self, indeg): + assert isinstance(indeg, bool), "`indeg` must be a boolean." + self.indeg = indeg + BaseRVS.__init__(self) + + def __call__(self, nsamples, random_state=42, dtype=numpy.float32): + gen = numpy.random.default_rng(random_state) + ra = 2 * numpy.pi * gen.random(nsamples, dtype=dtype) + dec = numpy.arcsin(2 * (gen.random(nsamples, dtype=dtype) - 0.5)) + if self.indeg: + ra = numpy.rad2deg(ra) + dec = numpy.rad2deg(dec) + return numpy.vstack([ra, dec]).T + + +############################################################################### +# RA wrapping # +############################################################################### + + +def wrapRA(ra, indeg): + """ + Wrap RA from :math:`[-180, 180)` to :math`[0, 360)` degrees if `indeg` or + equivalently in radians otherwise. + + Paramaters + ---------- + ra : 1-dimensional array + Right ascension. + indeg : bool + Whether the right ascension is in degrees. + + Returns + ------- + wrapped_ra : 1-dimensional array + """ + mask = ra < 0 + if numpy.sum(mask) == 0: + warn("No negative right ascension found.", UserWarning()) + ra[mask] += 360 if indeg else 2 * numpy.pi + return ra + + +############################################################################### +# Secondary assembly bias normalised marks # +############################################################################### + + +def normalised_marks(x, y, nbins): + """ + Calculate the normalised marks of `y` binned by `x`. + + Parameters + ---------- + x : 1-dimensional array + Binning variable. + y : 1-dimensional array + The variable to be marked. + nbins : int + Number of percentile bins. + + Returns + ------- + marks : 1-dimensional array + """ + assert x.ndim == y.ndim == 1 + if y.dtype not in [numpy.float32, numpy.float64]: + raise NotImplemented("Marks from integers are not supported.") + + bins = numpy.percentile(x, q=numpy.linspace(0, 100, nbins + 1)) + marks = numpy.full_like(y, numpy.nan) + for i in range(nbins): + m = (x >= bins[i]) & (x < bins[i + 1]) + # Calculate the normalised marks of this bin + _marks = numpy.full(numpy.sum(m), numpy.nan, dtype=marks.dtype) + for n, ind in enumerate(numpy.argsort(y[m])): + _marks[ind] = n + _marks /= numpy.nanmax(_marks) + marks[m] = _marks + + return marks diff --git a/csiborgtools/match/__init__.py b/csiborgtools/match/__init__.py index 4f7ac07..40c9e68 100644 --- a/csiborgtools/match/__init__.py +++ b/csiborgtools/match/__init__.py @@ -18,5 +18,3 @@ from .match import (RealisationsMatcher, cosine_similarity, # noqa calculate_overlap, calculate_overlap_indxs, # noqa dist_centmass, dist_percentile) # noqa from .num_density import (binned_counts, number_density) # noqa -from .knn import kNN_CDF -# from .correlation import (get_randoms_sphere, sphere_angular_tpcf) # noqa diff --git a/csiborgtools/read/summaries.py b/csiborgtools/read/summaries.py index 86e22d7..b2f37da 100644 --- a/csiborgtools/read/summaries.py +++ b/csiborgtools/read/summaries.py @@ -18,6 +18,7 @@ Tools for summarising various results. from os.path import (join, isfile) from glob import glob import numpy +from scipy.special import factorial import joblib from tqdm import tqdm @@ -184,55 +185,53 @@ class kNNCDFReader: """ Shortcut object to read in the kNN CDF data. """ - def read(self, files, ks, rmin=None, rmax=None, to_clip=True): + def read(self, run, folder, rmin=None, rmax=None, to_clip=True): """ - Read the kNN CDF data can be either the auto- or cross-correlation. + Read the auto- or cross-correlation kNN-CDF data. Infers the type from + the data files. Parameters ---------- - files : list of str - List of file paths to read in. - ks : list of int - kNN values to read in. + run : str + Run ID to read in. + folder : str + Path to the folder where the auto-correlation kNN-CDF is stored. rmin : float, optional Minimum separation. By default ignored. rmax : float, optional Maximum separation. By default ignored. to_clip : bool, optional - Whether to clip the auto-correlation CDF. Ignored if reading in the + Whether to clip the auto-correlation CDF. Ignored for cross-correlation. Returns ------- - rs : 1-dimensional array - Array of separations. - out : 4-dimensional array - Auto-correlation or cross-correlation kNN CDFs. The shape is - `(len(files), len(mass_thresholds), len(ks), neval)`. - mass_thresholds : 1-dimensional array - Array of mass thresholds. + rs : 1-dimensional array of shape `(neval, )` + Separations where the CDF is evaluated. + out : 3-dimensional array of shape `(len(files), len(ks), neval)` + Array of CDFs or cross-correlations. """ - data = joblib.load(files[0]) - if "cdf_0" in data.keys(): - isauto = True - kind = "cdf" - elif "corr_0" in data.keys(): - isauto = False - kind = "corr" - else: - raise ValueError("Unknown data format.") - rs = data["rs"] - mass_thresholds = data["mass_threshold"] - neval = data["{}_0".format(kind)].shape[1] - out = numpy.full((len(files), len(mass_thresholds), len(ks), neval), - numpy.nan, dtype=numpy.float32) + run += ".p" + files = [f for f in glob(join(folder, "*")) if run in f] + if len(files) == 0: + raise RuntimeError("No files found for run `{}`.".format(run[:-2])) - for i, file in enumerate(tqdm(files)): + for i, file in enumerate(files): data = joblib.load(file) - for j in range(len(mass_thresholds)): - out[i, j, ...] = data["{}_{}".format(kind, j)][ks, :] - if isauto and to_clip: - out[i, j, ...] = self.clipped_cdf(out[i, j, ...]) + if i == 0: # Initialise the array + if "corr" in data.keys(): + kind = "corr" + isauto = False + else: + kind = "cdf" + isauto = True + out = numpy.full((len(files), *data[kind].shape), numpy.nan, + dtype=numpy.float32) + rs = data["rs"] + out[i, ...] = data[kind] + + if isauto and to_clip: + out[i, ...] = self.clipped_cdf(out[i, ...]) # Apply separation cuts mask = (rs >= rmin if rmin is not None else rs > 0) @@ -240,7 +239,7 @@ class kNNCDFReader: rs = rs[mask] out = out[..., mask] - return rs, out, mass_thresholds + return rs, out @staticmethod def peaked_cdf(cdf, make_copy=True): @@ -295,37 +294,74 @@ class kNNCDFReader: return cdf @staticmethod - def prob_kvolume(cdfs, rs=None, normalise=False): - """ - Calculate the probability that a spherical volume contains :math:`k`= - objects from the kNN CDFs. + def prob_k(cdf): + r""" + Calculate the PDF that a spherical volume of radius :math:`r` contains + :math:`k` objects, i.e. :math:`P(k | V = 4 \pi r^3 / 3)`. Parameters ---------- - cdf : 4-dimensional array of shape `(nfiles, nmasses, nknn, nrs)` + cdf : 3-dimensional array of shape `(len(files), len(ks), len(rs))` Array of CDFs - normalise : bool, optional - Whether to normalise the probability to 1. Returns ------- - pk : 4-dimensional array of shape `(nfiles, nmasses, nknn - 1, nrs)` + pk : 3-dimensional array of shape `(len(files), len(ks)- 1, len(rs))` """ - out = numpy.full_like(cdfs[..., 1:, :], numpy.nan, dtype=numpy.float32) + out = numpy.full_like(cdf[..., 1:, :], numpy.nan, dtype=numpy.float32) + nks = cdf.shape[-2] + out[..., 0, :] = 1 - cdf[..., 0, :] - for k in range(cdfs.shape[-2] - 1): - out[..., k, :] = cdfs[..., k, :] - cdfs[..., k + 1, :] + for k in range(1, nks - 1): + out[..., k, :] = cdf[..., k - 1, :] - cdf[..., k, :] - if normalise: - assert rs is not None, "rs must be provided to normalise." - assert rs.ndim == 1 - - norm = numpy.nansum( - 0.5 * (out[..., 1:] + out[..., :-1]) * (rs[1:] - rs[:-1]), - axis=-1) - out /= norm.reshape(*norm.shape, 1) return out + def mean_prob_k(self, cdf): + """ + Calculate the mean PDF that a spherical volume of radius :math:`r` + contains :math:`k` objects, i.e. :math:`P(k | V = 4 \pi r^3 / 3)`, + averaged over the IC realisations. + + Parameters + ---------- + cdf : 3-dimensional array of shape `(len(files), len(ks), len(rs))` + Array of CDFs + Returns + ------- + out : 3-dimensional array of shape `(len(ks) - 1, len(rs), 2)` + Mean :math:`P(k | V = 4 \pi r^3 / 3) and its standard deviation, + stored along the last dimension, respectively. + """ + pk = self.prob_k(cdf) + return numpy.stack([numpy.mean(pk, axis=0), numpy.std(pk, axis=0)], + axis=-1) + + def poisson_prob_k(self, rs, k, ndensity): + """ + Calculate the analytical PDF that a spherical volume of + radius :math:`r` contains :math:`k` objects, i.e. + :math:`P(k | V = 4 \pi r^3 / 3)`, assuming a Poisson field (uniform + distribution of points). + + Parameters + ---------- + rs : 1-dimensional array + Array of separations. + k : int + Number of objects. + ndensity : float + Number density of objects. + + Returns + ------- + pk : 1-dimensional array + The PDF that a spherical volume of radius :math:`r` contains + :math:`k` objects. + """ + V = 4 * numpy.pi / 3 * rs**3 + return (ndensity * V)**k / factorial(k) * numpy.exp(-ndensity * V) + @staticmethod def cross_files(ic, folder): """ diff --git a/notebooks/knn.ipynb b/notebooks/knn.ipynb index ce76a00..7dd7ec5 100644 --- a/notebooks/knn.ipynb +++ b/notebooks/knn.ipynb @@ -6,8 +6,8 @@ "id": "5a38ed25", "metadata": { "ExecuteTime": { - "end_time": "2023-04-05T10:08:35.127670Z", - "start_time": "2023-04-05T10:08:29.875068Z" + "end_time": "2023-04-09T17:12:51.060443Z", + "start_time": "2023-04-09T17:12:47.288759Z" }, "scrolled": true }, @@ -35,6 +35,7 @@ " import sys\n", " sys.path.append(\"../\")\n", " import csiborgtools\n", + "import yaml\n", "\n", "\n", "%matplotlib notebook\n", @@ -44,68 +45,50 @@ }, { "cell_type": "code", - "execution_count": 135, - "id": "3bbf38d6", + "execution_count": 44, + "id": "3d39187a", "metadata": { "ExecuteTime": { - "end_time": "2023-04-05T13:41:02.225736Z", - "start_time": "2023-04-05T13:41:00.875529Z" + "end_time": "2023-04-09T19:51:15.144264Z", + "start_time": "2023-04-09T19:51:14.672337Z" } }, "outputs": [], "source": [ - "knnreader = csiborgtools.read.kNNCDFReader()" - ] - }, - { - "cell_type": "code", - "execution_count": 142, - "id": "9a79dde6", - "metadata": { - "ExecuteTime": { - "end_time": "2023-04-05T13:42:26.370674Z", - "start_time": "2023-04-05T13:42:22.862756Z" - } - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "100%|██████████| 101/101 [00:01<00:00, 76.37it/s]\n" - ] - } - ], - "source": [ - "files = glob(\"/mnt/extraspace/rstiskalek/csiborg/knn/auto/*\")\n", + "with open('../scripts/knn_auto.yml', 'r') as file:\n", + " config = yaml.safe_load(file)\n", "\n", - "ks = [0, 1, 2, 3, 4, 5, 6, 7]\n", - "rs, cdf, thresholds = knnreader.read(files, ks, rmin=0.01, rmax=100)" + "paths = csiborgtools.read.CSiBORGPaths()\n", + "knnreader = csiborgtools.read.kNNCDFReader()\n", + "\n", + "auto_folder = \"/mnt/extraspace/rstiskalek/csiborg/knn/auto/\"\n", + "cross_folder = \"/mnt/extraspace/rstiskalek/csiborg/knn/cross/\"\n" ] }, { "cell_type": "code", - "execution_count": 143, - "id": "88de6882", + "execution_count": 52, + "id": "202cb57b", "metadata": { "ExecuteTime": { - "end_time": "2023-04-05T13:42:26.943147Z", - "start_time": "2023-04-05T13:42:26.372382Z" + "end_time": "2023-04-09T19:53:58.915750Z", + "start_time": "2023-04-09T19:53:55.595006Z" } }, "outputs": [], "source": [ - "pk = knnreader.prob_kvolume(cdf, rs, True)" + "rs, corr = knnreader.read(\"mass001\", cross_folder, rmin=1, rmax=50)\n", + "rs, corr_rand = knnreader.read(\"mass001_random\", auto_folder, rmin=1, rmax=50)" ] }, { "cell_type": "code", - "execution_count": 147, - "id": "dbb11ba2", + "execution_count": 60, + "id": "9b3adce7", "metadata": { "ExecuteTime": { - "end_time": "2023-04-05T13:43:45.754506Z", - "start_time": "2023-04-05T13:43:38.953908Z" + "end_time": "2023-04-09T19:55:25.114871Z", + "start_time": "2023-04-09T19:55:24.983841Z" } }, "outputs": [ @@ -1077,7 +1060,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -1087,6 +1070,1772 @@ "output_type": "display_data" } ], + "source": [ + "cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n", + "plt.figure()\n", + "for k in range(64):\n", + " plt.plot(rs, corr[0, k, :], c=cols[0], lw=0.5)\n", + "# plt.plot(rs, np.abs(corr_rand[0, k, :]), c=cols[1], lw=0.5)\n", + "\n", + "plt.xscale(\"log\")\n", + "plt.axvline(2.65 / 0.705, c=\"red\", ls=\"--\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ec93a98", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "952cc374", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e38b7cf9", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "d46e176f", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-09T16:12:21.080837Z", + "start_time": "2023-04-09T16:12:19.137559Z" + } + }, + "outputs": [ + { + "ename": "RuntimeError", + "evalue": "No files found for run `mass001_spinloww`.", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn [57], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m rs, cdf \u001b[38;5;241m=\u001b[39m knnreader\u001b[38;5;241m.\u001b[39mread_auto(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmass001_spinloww\u001b[39m\u001b[38;5;124m\"\u001b[39m, folder, rmin\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m, rmax\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m50\u001b[39m)\n\u001b[1;32m 2\u001b[0m pk \u001b[38;5;241m=\u001b[39m knnreader\u001b[38;5;241m.\u001b[39mmean_prob_k(cdf)\n", + "File \u001b[0;32m~/csiborgtools/csiborgtools/read/summaries.py:215\u001b[0m, in \u001b[0;36mkNNCDFReader.read_auto\u001b[0;34m(self, run, folder, rmin, rmax, to_clip)\u001b[0m\n\u001b[1;32m 213\u001b[0m files \u001b[38;5;241m=\u001b[39m [f \u001b[38;5;28;01mfor\u001b[39;00m f \u001b[38;5;129;01min\u001b[39;00m glob(join(folder, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m*\u001b[39m\u001b[38;5;124m\"\u001b[39m)) \u001b[38;5;28;01mif\u001b[39;00m run \u001b[38;5;129;01min\u001b[39;00m f]\n\u001b[1;32m 214\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(files) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m--> 215\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mRuntimeError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNo files found for run `\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m`.\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(run[:\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m2\u001b[39m]))\n\u001b[1;32m 216\u001b[0m kind \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcorr\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcross\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01min\u001b[39;00m run \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcdf\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 218\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i, file \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28menumerate\u001b[39m(files):\n", + "\u001b[0;31mRuntimeError\u001b[0m: No files found for run `mass001_spinloww`." + ] + } + ], + "source": [ + "rs, cdf = knnreader.read_auto(\"mass001_spinloww\", folder, rmin=1, rmax=50)\n", + "pk = knnreader.mean_prob_k(cdf)\n", + "\n", + "\n", + "# rs, cdf = knnreader.read_auto(\"mass001_spinhigh\", folder, rmin=1, rmax=50)\n", + "# pk_perm = knnreader.mean_prob_k(cdf)" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "id": "ad36cab2", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-09T16:09:55.420866Z", + "start_time": "2023-04-09T16:09:55.391296Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(63, 5663, 2)" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pk.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "09015847", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-09T16:13:34.120967Z", + "start_time": "2023-04-09T16:13:31.837869Z" + } + }, + "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": [ + "cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n", + "\n", + "\n", + "plt.figure()\n", + "# plt.plot(rs, np.nansum(pk, axis=0)[:, 0], c=cols[k])\n", + "# plt.plot(rs, np.nansum(pk_perm, axis=0)[:, 0], c=cols[k], ls=\"--\")\n", + "for k in range(1, 2):\n", + " plt.plot(rs, pk[k, :, 0], c=cols[k])\n", + "# plt.plot(rs, pk_perm[k, :, 0], c=cols[k], ls=\"--\")\n", + "# plt.fill_between(rs, pk[k, :, 0] - pk[k, :, 1], pk[k, :, 0] + pk[k, :, 1])\n", + "\n", + "# plt.plot(rs, np.sum(pk, axis=0)[:, 0])\n", + "\n", + "# plt.xscale(\"log\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1279a8cb", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e0185ce", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "df973ed4", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T15:08:50.076207Z", + "start_time": "2023-04-08T15:08:46.683983Z" + }, + "scrolled": false + }, + "outputs": [], + "source": [ + "\n", + "\n", + "plt.figure()\n", + "\n", + "\n", + "k = 2\n", + "\n", + "rs, mu, std = mean_auto(k, \"mass_003_spinhigh\")\n", + "rs, mu_perm, std_perm = mean_auto(k, \"mass003_spinlow\")\n", + "z = mu / mu_perm\n", + "deltaz = z * np.sqrt((std / mu)**2 + (std_perm / mu_perm)**2)\n", + "plt.plot(rs, z, c=cols[0])\n", + "plt.fill_between(rs, z - deltaz, z + deltaz, color=cols[0], alpha=0.3)\n", + "\n", + "\n", + "# rs, mu, std = mean_auto(k, \"mass001_spinhigh\")\n", + "# rs, mu_perm, std_perm = mean_auto(k, \"mass001_spinhigh_perm\")\n", + "# z = mu / mu_perm\n", + "# deltaz = z * np.sqrt((std / mu)**2 + (std_perm / mu_perm)**2)\n", + "# plt.plot(rs, z, c=cols[1])\n", + "# plt.fill_between(rs, z - deltaz, z + deltaz, color=cols[1], alpha=0.3)\n", + "\n", + "\n", + "plt.axhline(1, c=\"black\", ls=\"--\")\n", + "# plt.fill_between(rs, mu - std, mu + std, color=cols[2], alpha=0.3)\n", + "\n", + "plt.xscale(\"log\")\n", + "# plt.yscale(\"log\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f500800", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T14:11:02.568517Z", + "start_time": "2023-04-08T14:11:00.406477Z" + } + }, + "outputs": [], + "source": [ + "# rs, corr = knnreader.read_auto(\"mass001_spinmedian_cross\", folder, rmin=0.5)\n", + "# rs, corr_perm = knnreader.read_auto(\"mass001_spinmedian_cross_perm\", folder, rmin=0.5)\n", + "\n", + "# rs, cdf_low = knnreader.read_auto(\"mass001_spinmedian_cross_perm\", folder, rmin=0.5)\n", + "rs, cdf_high = knnreader.read_auto(\"mass001_spinmedian_cross_perm\", folder, rmin=0.5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6565101a", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T14:04:17.514165Z", + "start_time": "2023-04-08T14:04:15.555281Z" + } + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84729726", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T15:04:01.025738Z", + "start_time": "2023-04-08T15:03:56.792105Z" + }, + "scrolled": false + }, + "outputs": [], + "source": [ + "plt.figure()\n", + "rs, mu, std = mean_auto(0, \"mass001_spinlow\")\n", + "plt.plot(rs, mu, c=cols[0])\n", + "plt.fill_between(rs, mu - std, mu + std, color=cols[0], alpha=0.5)\n", + "\n", + "rs, mu, std = mean_auto(0, \"mass001_spinlow_perm\")\n", + "plt.plot(rs, mu, c=cols[1])\n", + "plt.fill_between(rs, mu - std, mu + std, color=cols[1], alpha=0.5)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f85f943", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T14:11:05.066794Z", + "start_time": "2023-04-08T14:11:04.362662Z" + } + }, + "outputs": [], + "source": [ + "prk_low = knnreader.prk(rs, cdf_low)\n", + "prk_high = knnreader.prk(rs, cdf_high)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c862e8b", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T14:11:07.074775Z", + "start_time": "2023-04-08T14:11:05.148559Z" + } + }, + "outputs": [], + "source": [ + "k = 0\n", + "\n", + "plt.figure()\n", + "\n", + "for i in range(101):\n", + " plt.plot(rs, prk_low[i, k, :], lw=0.1, c=cols[0])\n", + " plt.plot(rs, prk_high[i, k, :], lw=0.1, c=cols[1])\n", + "\n", + " \n", + "rs, mu, std = mean_auto(0, \"mass001\")\n", + "plt.plot(rs, mu, c=cols[2])\n", + "plt.fill_between(rs, mu - std, mu + std, color=cols[2], alpha=0.5)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6682fa88", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T14:59:29.041395Z", + "start_time": "2023-04-08T14:59:24.875903Z" + } + }, + "outputs": [], + "source": [ + "rs, corr = knnreader.read_auto(\"mass001_spinmedian_cross\", folder, rmin=0.5)\n", + "rs, corrperm = knnreader.read_auto(\"mass001_spinmedian_cross_perm\", folder, rmin=0.5)\n", + "\n", + "# rs, corr_low = knnreader.read_auto(\"mass001_spinlow_cross_perm\", folder, rmin=0.5)\n", + "# rs, corr_high = knnreader.read_auto(\"mass001_spinhigh_cross_perm\", folder, rmin=0.5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "298f28b6", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T14:59:47.556944Z", + "start_time": "2023-04-08T14:59:47.103398Z" + } + }, + "outputs": [], + "source": [ + "k = 2\n", + "plt.figure()\n", + "for i in range(101):\n", + " plt.plot(rs, corr[i, k, :], lw=0.1, c=cols[0])\n", + " plt.plot(rs, corrperm[i, k, :], lw=0.1, c=cols[1])\n", + "# plt.plot(rs, corr[i, k, :] - corrperm[i, k, :], lw=0.1, c=cols[0])\n", + "# plt.plot(rs, , lw=0.1, c=cols[1])\n", + " \n", + "plt.xscale(\"log\")\n", + "plt.yscale(\"log\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a28d0290", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T13:59:25.401233Z", + "start_time": "2023-04-08T13:59:25.157683Z" + } + }, + "outputs": [], + "source": [ + "k = 0\n", + "\n", + "plt.figure()\n", + "for i in range(101):\n", + " plt.plot(rs, corr[i, k, :], c=cols[0], lw=0.1)\n", + " \n", + " \n", + "for i in range(101):\n", + " plt.plot(rs, corr_perm[i, k, :], c=cols[1], lw=0.1)\n", + "\n", + "plt.xscale(\"log\")\n", + "plt.yscale(\"log\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e86ce611", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-08T09:30:44.330451Z", + "start_time": "2023-04-08T09:30:44.295713Z" + } + }, + "outputs": [], + "source": [ + "mean_prk.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe9a0ef3", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "de17207d", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4cbad0c", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe11b032", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "491be45e", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f19accd3", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T16:51:43.224223Z", + "start_time": "2023-04-07T16:51:35.597574Z" + } + }, + "outputs": [], + "source": [ + "cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n", + "\n", + "\n", + "plt.figure()\n", + "\n", + "rs, cdf = knnreader.read_auto(\"mass001_spinlow\", folder, rmin=0.5)\n", + "pk = knnreader.prob_kvolume(cdf, rs, True)\n", + "for i in range(101):\n", + " plt.plot(rs, pk[i, 5, :], lw=0.1, c=cols[0])\n", + " \n", + " \n", + "rs, cdf = knnreader.read_auto(\"mass001_spinhigh\", folder, rmin=0.5)\n", + "pk = knnreader.prob_kvolume(cdf, rs, True)\n", + "for i in range(101):\n", + " plt.plot(rs, pk[i, 5, :], lw=0.1, c=cols[1]) \n", + " \n", + " \n", + "rs, cdf = knnreader.read_auto(\"mass001_spinhigh_perm\", folder, rmin=0.5)\n", + "pk = knnreader.prob_kvolume(cdf, rs, True)\n", + "for i in range(101):\n", + " plt.plot(rs, pk[i, 5, :], lw=0.1, c=cols[2])\n", + "\n", + "# plt.xscale(\"log\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f735b44", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a273df53", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "900fd4f8", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01974708", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T12:20:08.521290Z", + "start_time": "2023-04-07T12:20:05.258954Z" + } + }, + "outputs": [], + "source": [ + "cat = csiborgtools.read.HaloCatalogue(7444, paths, min_mass=1e12, max_dist=155/0.705)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47c494f6", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86f02695", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T11:51:11.124583Z", + "start_time": "2023-04-07T11:51:11.094094Z" + } + }, + "outputs": [], + "source": [ + "x = \"\"\n", + "for key in auto_config.keys():\n", + " x += key + \" \"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "62d9e837", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T11:51:12.533128Z", + "start_time": "2023-04-07T11:51:12.499981Z" + } + }, + "outputs": [], + "source": [ + "x" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db6fc6e4", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T10:28:20.801576Z", + "start_time": "2023-04-07T10:28:20.766160Z" + } + }, + "outputs": [], + "source": [ + "auto_config" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8fdcfb01", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T10:21:14.858309Z", + "start_time": "2023-04-07T10:21:14.826448Z" + } + }, + "outputs": [], + "source": [ + "auto_folder = \"/mnt/extraspace/rstiskalek/csiborg/knn/auto/\"\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b5f37af", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T12:20:38.019809Z", + "start_time": "2023-04-07T12:20:37.987281Z" + } + }, + "outputs": [], + "source": [ + "np.log10(cat[\"totpartmass\"].min())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5b2df40", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T11:13:18.959385Z", + "start_time": "2023-04-07T11:13:13.072536Z" + } + }, + "outputs": [], + "source": [ + "cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n", + "\n", + "\n", + "plt.figure()\n", + "rs, cdf = knnreader.read_auto(\"004\", auto_folder)\n", + "pk = knnreader.prob_kvolume(cdf, rs, normalise=True)\n", + "for i in range(101):\n", + " plt.plot(rs, pk[i, 0, :], c=cols[0], lw=0.1)\n", + "\n", + "\n", + "rs, cdf = knnreader.read_auto(\"005\", auto_folder)\n", + "pk = knnreader.prob_kvolume(cdf, rs, normalise=True)\n", + "for i in range(101):\n", + " plt.plot(rs, pk[i, 0, :], c=cols[1], lw=0.1)\n", + "\n", + " \n", + "rs, cdf = knnreader.read_auto(\"001\", auto_folder)\n", + "pk = knnreader.prob_kvolume(cdf, rs, normalise=True)\n", + "for i in range(101):\n", + " plt.plot(rs, pk[i, 0, :], c=cols[2], lw=0.1)\n", + "\n", + "# plt.xscale(\"log\")\n", + "# plt.yscale(\"log\")\n", + "\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "85c65eef", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20ecf551", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ac70147", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T09:47:03.873611Z", + "start_time": "2023-04-07T09:47:01.794363Z" + } + }, + "outputs": [], + "source": [ + "cat = csiborgtools.read.HaloCatalogue(7444, paths)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b77b1377", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T07:49:23.539504Z", + "start_time": "2023-04-07T07:45:50.514216Z" + } + }, + "outputs": [], + "source": [ + "from tqdm import trange\n", + "x = np.full((len(ics), 3), np.nan)\n", + "for i in trange(len(ics)):\n", + " cat = csiborgtools.read.HaloCatalogue(ics[i], paths, max_dist=155 / 0.705)\n", + " for j, th in enumerate([1e12, 1e13, 1e14]):\n", + " mask = cat[\"totpartmass\"] > th\n", + " x[i, j] = np.nanmedian(cat[\"lambda200c\"][mask])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e55d821d", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T07:51:55.598449Z", + "start_time": "2023-04-07T07:51:55.567861Z" + } + }, + "outputs": [], + "source": [ + "np.mean(x[:, 2]), np.std(x[:, 2])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3bae7373", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-07T11:12:14.971055Z", + "start_time": "2023-04-07T11:12:14.938117Z" + } + }, + "outputs": [], + "source": [ + "cdf.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "127cb78e", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8368d476", + "metadata": {}, + "outputs": [], + "source": [ + "np.nanmedian()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6b5a68c", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-06T17:59:09.028319Z", + "start_time": "2023-04-06T17:59:08.953783Z" + } + }, + "outputs": [], + "source": [ + "plt.figure()\n", + "\n", + "plt.scatter(cat[\"m200\"], cat[\"lambda200c\"], s=1)\n", + "\n", + "plt.xscale(\"log\")\n", + "plt.yscale(\"log\")\n", + "\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2eb1c2e4", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a79dde6", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-05T13:42:26.370674Z", + "start_time": "2023-04-05T13:42:22.862756Z" + } + }, + "outputs": [], + "source": [ + "files = glob(\"/mnt/extraspace/rstiskalek/csiborg/knn/auto/*\")\n", + "\n", + "ks = [0, 1, 2, 3, 4, 5, 6, 7]\n", + "rs, cdf, thresholds = knnreader.read(files, ks, rmin=0.01, rmax=100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88de6882", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-05T13:42:26.943147Z", + "start_time": "2023-04-05T13:42:26.372382Z" + } + }, + "outputs": [], + "source": [ + "pk = knnreader.prob_kvolume(cdf, rs, True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dbb11ba2", + "metadata": { + "ExecuteTime": { + "end_time": "2023-04-05T13:43:45.754506Z", + "start_time": "2023-04-05T13:43:38.953908Z" + } + }, + "outputs": [], "source": [ "cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n", "\n", diff --git a/scripts/run_crosspk.py b/scripts/crosspk.py similarity index 100% rename from scripts/run_crosspk.py rename to scripts/crosspk.py diff --git a/scripts/run_fieldprop.py b/scripts/fieldprop.py similarity index 100% rename from scripts/run_fieldprop.py rename to scripts/fieldprop.py diff --git a/scripts/run_fit_halos.py b/scripts/fit_halos.py similarity index 100% rename from scripts/run_fit_halos.py rename to scripts/fit_halos.py diff --git a/scripts/run_initmatch.py b/scripts/initmatch.py similarity index 100% rename from scripts/run_initmatch.py rename to scripts/initmatch.py diff --git a/scripts/knn_auto.py b/scripts/knn_auto.py new file mode 100644 index 0000000..2dc59bf --- /dev/null +++ b/scripts/knn_auto.py @@ -0,0 +1,182 @@ +# 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. +"""A script to calculate the KNN-CDF for a set of CSiBORG halo catalogues.""" +from os.path import join +from warnings import warn +from argparse import ArgumentParser +from copy import deepcopy +from datetime import datetime +from mpi4py import MPI +from TaskmasterMPI import master_process, worker_process +import numpy +from sklearn.neighbors import NearestNeighbors +import joblib +import yaml +try: + import csiborgtools +except ModuleNotFoundError: + import sys + sys.path.append("../") + import csiborgtools + + +############################################################################### +# MPI and arguments # +############################################################################### +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +nproc = comm.Get_size() + +parser = ArgumentParser() +parser.add_argument("--runs", type=str, nargs="+") +args = parser.parse_args() +with open('../scripts/knn_auto.yml', 'r') as file: + config = yaml.safe_load(file) + +Rmax = 155 / 0.705 # Mpc (h = 0.705) high resolution region radius +totvol = 4 * numpy.pi * Rmax**3 / 3 +minmass = 1e12 +ics = [7444, 7468, 7492, 7516, 7540, 7564, 7588, 7612, 7636, 7660, 7684, + 7708, 7732, 7756, 7780, 7804, 7828, 7852, 7876, 7900, 7924, 7948, + 7972, 7996, 8020, 8044, 8068, 8092, 8116, 8140, 8164, 8188, 8212, + 8236, 8260, 8284, 8308, 8332, 8356, 8380, 8404, 8428, 8452, 8476, + 8500, 8524, 8548, 8572, 8596, 8620, 8644, 8668, 8692, 8716, 8740, + 8764, 8788, 8812, 8836, 8860, 8884, 8908, 8932, 8956, 8980, 9004, + 9028, 9052, 9076, 9100, 9124, 9148, 9172, 9196, 9220, 9244, 9268, + 9292, 9316, 9340, 9364, 9388, 9412, 9436, 9460, 9484, 9508, 9532, + 9556, 9580, 9604, 9628, 9652, 9676, 9700, 9724, 9748, 9772, 9796, + 9820, 9844] +dumpdir = "/mnt/extraspace/rstiskalek/csiborg/knn" +fout = join(dumpdir, "auto", "knncdf_{}_{}.p") +paths = csiborgtools.read.CSiBORGPaths() +knncdf = csiborgtools.clustering.kNN_CDF() + +############################################################################### +# Analysis # +############################################################################### + +def read_single(selection, cat): + """Positions for single catalogue auto-correlation.""" + mmask = numpy.ones(len(cat), dtype=bool) + pos = cat.positions(False) + # Primary selection + psel = selection["primary"] + pmin, pmax = psel.get("min", None), psel.get("max", None) + if pmin is not None: + mmask &= (cat[psel["name"]] >= pmin) + if pmax is not None: + mmask &= (cat[psel["name"]] < pmax) + pos = pos[mmask, ...] + + # Secondary selection + if "secondary" not in selection: + return pos + smask = numpy.ones(pos.shape[0], dtype=bool) + ssel = selection["secondary"] + smin, smax = ssel.get("min", None), ssel.get("max", None) + prop = cat[ssel["name"]][mmask] + if ssel.get("toperm", False): + prop = numpy.random.permutation(prop) + if ssel.get("marked", True): + x = cat[psel["name"]][mmask] + prop = csiborgtools.clustering.normalised_marks( + x, prop, nbins=config["nbins_marks"]) + + if smin is not None: + smask &= (prop >= smin) + if smax is not None: + smask &= (prop < smax) + + return pos[smask, ...] + +def do_auto(run, cat, ic): + """Calculate the kNN-CDF single catalgoue autocorrelation.""" + _config = config.get(run, None) + if _config is None: + warn("No configuration for run {}.".format(run)) + return + + rvs_gen = csiborgtools.clustering.RVSinsphere(Rmax) + pos = read_single(_config, cat) + knn = NearestNeighbors() + knn.fit(pos) + rs, cdf = knncdf( + knn, rvs_gen=rvs_gen, nneighbours=config["nneighbours"], + rmin=config["rmin"], rmax=config["rmax"], + nsamples=int(config["nsamples"]), neval=int(config["neval"]), + batch_size=int(config["batch_size"]), random_state=config["seed"]) + + joblib.dump({"rs": rs, "cdf": cdf, "ndensity": pos.shape[0] / totvol}, + fout.format(str(ic).zfill(5), run)) + +def do_cross_rand(run, cat, ic): + """Calculate the kNN-CDF cross catalogue random correlation.""" + _config = config.get(run, None) + if _config is None: + warn("No configuration for run {}.".format(run)) + return + + rvs_gen = csiborgtools.clustering.RVSinsphere(Rmax) + knn1, knn2 = NearestNeighbors(), NearestNeighbors() + + pos1 = read_single(_config, cat) + knn1.fit(pos1) + + pos2 = rvs_gen(pos1.shape[0]) + knn2.fit(pos2) + + rs, cdf0, cdf1, joint_cdf = knncdf.joint( + knn1, knn2, rvs_gen=rvs_gen, nneighbours=int(config["nneighbours"]), + rmin=config["rmin"], rmax=config["rmax"], + nsamples=int(config["nsamples"]), neval=int(config["neval"]), + batch_size=int(config["batch_size"]), random_state=config["seed"]) + corr = knncdf.joint_to_corr(cdf0, cdf1, joint_cdf) + + joblib.dump({"rs": rs, "corr": corr}, fout.format(str(ic).zfill(5), run)) + + + +def do_runs(ic): + cat = csiborgtools.read.HaloCatalogue(ic, paths, max_dist=Rmax, + min_mass=minmass) + for run in args.runs: + if "random" in run: + do_cross_rand(run, cat, ic) + else: + do_auto(run, cat, ic) + + +############################################################################### +# MPI task delegation # +############################################################################### + + +if nproc > 1: + if rank == 0: + tasks = deepcopy(ics) + master_process(tasks, comm, verbose=True) + else: + worker_process(do_runs, comm, verbose=False) +else: + tasks = deepcopy(ics) + for task in tasks: + print("{}: completing task `{}`.".format(datetime.now(), task)) + do_runs(task) +comm.Barrier() + + +if rank == 0: + print("{}: all finished.".format(datetime.now())) +quit() # Force quit the script diff --git a/scripts/knn_auto.yml b/scripts/knn_auto.yml new file mode 100644 index 0000000..ae02db9 --- /dev/null +++ b/scripts/knn_auto.yml @@ -0,0 +1,144 @@ +rmin: 0.1 +rmax: 100 +nneighbours: 64 +nsamples: 1.e+7 +batch_size: 1.e+6 +neval: 10000 +seed: 42 +nbins_marks: 10 + + +################################################################################ +# totpartmass # +################################################################################ + + +"mass001": + primary: + name: totpartmass + min: 1.e+12 + max: 1.e+13 + +"mass002": + primary: + name: totpartmass + min: 1.e+13 + max: 1.e+14 + +"mass003": + primary: + name: totpartmass + min: 1.e+14 + + +################################################################################ +# totpartmass + lambda200c # +################################################################################ + + +"mass001_spinlow": + primary: + name: totpartmass + min: 1.e+12 + max: 1.e+13 + secondary: + name: lambda200c + toperm: false + marked: false + max: 0.5 + +"mass001_spinhigh": + primary: + name: totpartmass + min: 1.e+12 + max: 1.e+13 + secondary: + name: lambda200c + toperm: false + marked: true + min: 0.5 + +"mass001_spinmedian_perm": + primary: + name: totpartmass + min: 1.e+12 + max: 1.e+13 + secondary: + name: lambda200c + toperm: true + marked : true + min: 0.5 + +"mass002_spinlow": + primary: + name: totpartmass + min: 1.e+13 + max: 1.e+14 + secondary: + name: lambda200c + toperm: false + marked: false + max: 0.5 + +"mass002_spinhigh": + primary: + name: totpartmass + min: 1.e+13 + max: 1.e+14 + secondary: + name: lambda200c + toperm: false + marked: true + min: 0.5 + +"mass002_spinmedian_perm": + primary: + name: totpartmass + min: 1.e+13 + max: 1.e+14 + secondary: + name: lambda200c + toperm: true + marked : true + min: 0.5 + +"mass003_spinlow": + primary: + name: totpartmass + min: 1.e+14 + secondary: + name: lambda200c + toperm: false + marked: false + max: 0.5 + +"mass003_spinhigh": + primary: + name: totpartmass + min: 1.e+14 + secondary: + name: lambda200c + toperm: false + marked: true + min: 0.5 + +"mass003_spinmedian_perm": + primary: + name: totpartmass + min: 1.e+14 + secondary: + name: lambda200c + toperm: true + marked : true + min: 0.5 + + +################################################################################ +# Cross with random # +################################################################################ + +"mass001_random": + primary: + name: totpartmass + min: 1.e+12 + max: 1.e+13 \ No newline at end of file diff --git a/scripts/run_knn.py b/scripts/knn_cross.py similarity index 55% rename from scripts/run_knn.py rename to scripts/knn_cross.py index d591eb3..fbccfed 100644 --- a/scripts/run_knn.py +++ b/scripts/knn_cross.py @@ -13,6 +13,7 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. """A script to calculate the KNN-CDF for a set of CSiBORG halo catalogues.""" +from warnings import warn from os.path import join from argparse import ArgumentParser from copy import deepcopy @@ -20,8 +21,10 @@ from datetime import datetime from itertools import combinations from mpi4py import MPI from TaskmasterMPI import master_process, worker_process +import numpy from sklearn.neighbors import NearestNeighbors import joblib +import yaml try: import csiborgtools except ModuleNotFoundError: @@ -38,17 +41,13 @@ rank = comm.Get_rank() nproc = comm.Get_size() parser = ArgumentParser() -parser.add_argument("--rmin", type=float) -parser.add_argument("--rmax", type=float) -parser.add_argument("--nneighbours", type=int) -parser.add_argument("--nsamples", type=int) -parser.add_argument("--neval", type=int) -parser.add_argument("--batch_size", type=int) -parser.add_argument("--seed", type=int, default=42) +parser.add_argument("--runs", type=str, nargs="+") args = parser.parse_args() +with open('../scripts/knn_cross.yml', 'r') as file: + config = yaml.safe_load(file) -Rmax = 155 / 0.705 # Mpc/h high resolution region radius -mass_threshold = [1e12, 1e13, 1e14] # Msun +Rmax = 155 / 0.705 # Mpc (h = 0.705) high resolution region radius +minmass = 1e12 ics = [7444, 7468, 7492, 7516, 7540, 7564, 7588, 7612, 7636, 7660, 7684, 7708, 7732, 7756, 7780, 7804, 7828, 7852, 7876, 7900, 7924, 7948, 7972, 7996, 8020, 8044, 8068, 8092, 8116, 8140, 8164, 8188, 8212, @@ -59,80 +58,58 @@ ics = [7444, 7468, 7492, 7516, 7540, 7564, 7588, 7612, 7636, 7660, 7684, 9292, 9316, 9340, 9364, 9388, 9412, 9436, 9460, 9484, 9508, 9532, 9556, 9580, 9604, 9628, 9652, 9676, 9700, 9724, 9748, 9772, 9796, 9820, 9844] -dumpdir = "/mnt/extraspace/rstiskalek/csiborg/knn" -fout_auto = join(dumpdir, "auto", "knncdf_{}.p") -fout_cross = join(dumpdir, "cross", "knncdf_{}_{}.p") paths = csiborgtools.read.CSiBORGPaths() - +dumpdir = "/mnt/extraspace/rstiskalek/csiborg/knn" +fout = join(dumpdir, "cross", "knncdf_{}_{}_{}.p") +knncdf = csiborgtools.clustering.kNN_CDF() ############################################################################### # Analysis # ############################################################################### -knncdf = csiborgtools.match.kNN_CDF() +def read_single(selection, cat): + mmask = numpy.ones(len(cat), dtype=bool) + pos = cat.positions(False) + # Primary selection + psel = selection["primary"] + pmin, pmax = psel.get("min", None), psel.get("max", None) + if pmin is not None: + mmask &= (cat[psel["name"]] >= pmin) + if pmax is not None: + mmask &= (cat[psel["name"]] < pmax) + return pos[mmask, ...] -def do_auto(ic): - out = {} - cat = csiborgtools.read.HaloCatalogue(ic, paths, max_dist=Rmax) +def do_cross(run, ics): + _config = config.get(run, None) + if _config is None: + warn("No configuration for run {}.".format(run)) + return + rvs_gen = csiborgtools.clustering.RVSinsphere(Rmax) + knn1, knn2 = NearestNeighbors(), NearestNeighbors() - for i, mmin in enumerate(mass_threshold): - knn = NearestNeighbors() - knn.fit(cat.positions(False)[cat["totpartmass"] > mmin, ...]) - - rs, cdf = knncdf(knn, nneighbours=args.nneighbours, Rmax=Rmax, - rmin=args.rmin, rmax=args.rmax, nsamples=args.nsamples, - neval=args.neval, batch_size=args.batch_size, - random_state=args.seed, verbose=False) - out.update({"cdf_{}".format(i): cdf}) - - out.update({"rs": rs, "mass_threshold": mass_threshold}) - joblib.dump(out, fout_auto.format(ic)) - - -def do_cross(ics): - out = {} cat1 = csiborgtools.read.HaloCatalogue(ics[0], paths, max_dist=Rmax) + pos1 = read_single(_config, cat1) + knn1.fit(pos1) + cat2 = csiborgtools.read.HaloCatalogue(ics[1], paths, max_dist=Rmax) + pos2 = read_single(_config, cat2) + knn2.fit(pos2) - for i, mmin in enumerate(mass_threshold): - knn1 = NearestNeighbors() - knn1.fit(cat1.positions()[cat1["totpartmass"] > mmin, ...]) + rs, cdf0, cdf1, joint_cdf = knncdf.joint( + knn1, knn2, rvs_gen=rvs_gen, nneighbours=int(config["nneighbours"]), + rmin=config["rmin"], rmax=config["rmax"], + nsamples=int(config["nsamples"]), neval=int(config["neval"]), + batch_size=int(config["batch_size"]), random_state=config["seed"]) - knn2 = NearestNeighbors() - knn2.fit(cat2.positions()[cat2["totpartmass"] > mmin, ...]) + corr = knncdf.joint_to_corr(cdf0, cdf1, joint_cdf) - rs, cdf0, cdf1, joint_cdf = knncdf.joint( - knn1, knn2, nneighbours=args.nneighbours, Rmax=Rmax, - rmin=args.rmin, rmax=args.rmax, nsamples=args.nsamples, - neval=args.neval, batch_size=args.batch_size, - random_state=args.seed) + joblib.dump({"rs": rs, "corr": corr}, + fout.format(str(ics[0]).zfill(5), str(ics[1]).zfill(5), run)) - corr = knncdf.joint_to_corr(cdf0, cdf1, joint_cdf) - - out.update({"corr_{}".format(i): corr}) - - out.update({"rs": rs, "mass_threshold": mass_threshold}) - joblib.dump(out, fout_cross.format(*ics)) - - - -############################################################################### -# Autocorrelation calculation # -############################################################################### - - -if nproc > 1: - if rank == 0: - tasks = deepcopy(ics) - master_process(tasks, comm, verbose=True) - else: - worker_process(do_auto, comm, verbose=False) -else: - tasks = deepcopy(ics) - for task in tasks: - print("{}: completing task `{}`.".format(datetime.now(), task)) - do_auto(task) -comm.Barrier() +def do_runs(ics): + print(ics) + for run in args.runs: + do_cross(run, ics) ############################################################################### @@ -145,12 +122,12 @@ if nproc > 1: tasks = list(combinations(ics, 2)) master_process(tasks, comm, verbose=True) else: - worker_process(do_cross, comm, verbose=False) + worker_process(do_runs, comm, verbose=False) else: - tasks = deepcopy(ics) + tasks = list(combinations(ics, 2)) for task in tasks: print("{}: completing task `{}`.".format(datetime.now(), task)) - do_cross(task) + do_runs(task) comm.Barrier() diff --git a/scripts/knn_cross.yml b/scripts/knn_cross.yml new file mode 100644 index 0000000..53e214a --- /dev/null +++ b/scripts/knn_cross.yml @@ -0,0 +1,29 @@ +rmin: 0.1 +rmax: 100 +nneighbours: 64 +nsamples: 1.e+7 +batch_size: 1.e+6 +neval: 10000 +seed: 42 + + +################################################################################ +# totpartmass # +################################################################################ + +"mass001": + primary: + name: totpartmass + min: 1.e+12 + max: 1.e+13 + +"mass002": + primary: + name: totpartmass + min: 1.e+13 + max: 1.e+14 + +"mass003": + primary: + name: totpartmass + min: 1.e+14 diff --git a/scripts/python.sh b/scripts/python.sh deleted file mode 100644 index 45328c4..0000000 --- a/scripts/python.sh +++ /dev/null @@ -1,46 +0,0 @@ -#!/bin/bash -l -echo ========================================================= -echo Job submitted date = Fri Mar 31 16:17:57 BST 2023 -date_start=`date +%s` -echo $SLURM_JOB_NUM_NODES nodes \( $SMP processes per node \) -echo $SLURM_JOB_NUM_NODES hosts used: $SLURM_JOB_NODELIST -echo Job output begins -echo ----------------- -echo -#hostname - -# Need to set the max locked memory very high otherwise IB can't allocate enough and fails with "UCX ERROR Failed to allocate memory pool chunk: Input/output error" -ulimit -l unlimited - -# To allow mvapich to run ok -export MV2_SMP_USE_CMA=0 - -#which mpirun -export OMP_NUM_THEADS=1 - /usr/local/shared/slurm/bin/srun -u -n 5 --mpi=pmi2 --mem-per-cpu=7168 nice -n 10 /mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python run_knn.py --rmin 0.05 --rmax 50 --nsamples 100000 --neval 10000 -# If we've been checkpointed -#if [ -n "${DMTCP_CHECKPOINT_DIR}" ]; then - if [ -d "${DMTCP_CHECKPOINT_DIR}" ]; then -# echo -n "Job was checkpointed at " -# date -# echo - sleep 1 -# fi - echo -n -else - echo --------------- - echo Job output ends - date_end=`date +%s` - seconds=$((date_end-date_start)) - minutes=$((seconds/60)) - seconds=$((seconds-60*minutes)) - hours=$((minutes/60)) - minutes=$((minutes-60*hours)) - echo ========================================================= - echo PBS job: finished date = `date` - echo Total run time : $hours Hours $minutes Minutes $seconds Seconds - echo ========================================================= -fi -if [ ${SLURM_NTASKS} -eq 1 ]; then - rm -f $fname -fi diff --git a/scripts/run_crosspk.sh b/scripts/run_crosspk.sh deleted file mode 100644 index 374cb88..0000000 --- a/scripts/run_crosspk.sh +++ /dev/null @@ -1,14 +0,0 @@ -nthreads=20 -memory=40 -queue="berg" -env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python" -file="run_crosspk.py" -grid=1024 -halfwidth=0.13 - -cm="addqueue -q $queue -n $nthreads -m $memory $env $file --grid $grid --halfwidth $halfwidth" - -echo "Submitting:" -echo $cm -echo -$cm diff --git a/scripts/run_fieldprop.sh b/scripts/run_fieldprop.sh deleted file mode 100644 index 76e8e17..0000000 --- a/scripts/run_fieldprop.sh +++ /dev/null @@ -1,14 +0,0 @@ -nthreads=10 -memory=32 -queue="berg" -env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python" -file="run_fieldprop.py" -# grid=1024 -# halfwidth=0.1 - -cm="addqueue -q $queue -n $nthreads -m $memory $env $file" - -echo "Submitting:" -echo $cm -echo -$cm diff --git a/scripts/run_fit_halos.sh b/scripts/run_fit_halos.sh deleted file mode 100644 index 5bc4f8b..0000000 --- a/scripts/run_fit_halos.sh +++ /dev/null @@ -1,12 +0,0 @@ -nthreads=100 -memory=3 -queue="berg" -env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python" -file="run_fit_halos.py" - -cm="addqueue -q $queue -n $nthreads -m $memory $env $file" - -echo "Submitting:" -echo $cm -echo -$cm diff --git a/scripts/run_initmatch.sh b/scripts/run_initmatch.sh deleted file mode 100644 index 3de6233..0000000 --- a/scripts/run_initmatch.sh +++ /dev/null @@ -1,14 +0,0 @@ -nthreads=15 # There isn't too much benefit going to too many CPUs... -memory=32 -queue="berg" -env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python" -file="run_initmatch.py" - -dump_clumps="false" - -cm="addqueue -q $queue -n $nthreads -m $memory $env $file --dump_clumps $dump_clumps" - -echo "Submitting:" -echo $cm -echo -$cm diff --git a/scripts/run_knn.sh b/scripts/run_knn.sh deleted file mode 100644 index a0c84a3..0000000 --- a/scripts/run_knn.sh +++ /dev/null @@ -1,23 +0,0 @@ -nthreads=151 -memory=4 -queue="cmb" -env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python" -file="run_knn.py" - -rmin=0.01 -rmax=100 -nneighbours=8 -nsamples=100000000 -batch_size=1000000 -neval=10000 - -pythoncm="$env $file --rmin $rmin --rmax $rmax --nneighbours $nneighbours --nsamples $nsamples --batch_size $batch_size --neval $neval" - -# echo $pythoncm -# $pythoncm - -cm="addqueue -q $queue -n $nthreads -m $memory $pythoncm" -echo "Submitting:" -echo $cm -echo -$cm diff --git a/scripts/run_singlematch.sh b/scripts/run_singlematch.sh deleted file mode 100755 index 58cc8d7..0000000 --- a/scripts/run_singlematch.sh +++ /dev/null @@ -1,36 +0,0 @@ -#!/bin/bash -# nthreads=1 -memory=16 -queue="berg" -env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python" -file="run_singlematch.py" - -nmult=1. -sigma=1. - -sims=(7468 7588 8020 8452 8836) -nsims=${#sims[@]} - -for i in $(seq 0 $((nsims-1))); do -for j in $(seq 0 $((nsims-1))); do -if [ $i -eq $j ]; then - continue -elif [ $i -gt $j ]; then - continue -else - : -fi - -nsim0=${sims[$i]} -nsimx=${sims[$j]} - -pythoncm="$env $file --nsim0 $nsim0 --nsimx $nsimx --nmult $nmult --sigma $sigma" - -cm="addqueue -q $queue -n 1x1 -m $memory $pythoncm" -echo "Submitting:" -echo $cm -echo -$cm -sleep 0.05 - -done; done diff --git a/scripts/run_split_halos.sh b/scripts/run_split_halos.sh deleted file mode 100644 index 84e93af..0000000 --- a/scripts/run_split_halos.sh +++ /dev/null @@ -1,12 +0,0 @@ -nthreads=1 -memory=30 -queue="cmb" -env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python" -file="run_split_halos.py" - -cm="addqueue -q $queue -n $nthreads -m $memory $env $file" - -echo "Submitting:" -echo $cm -echo -$cm diff --git a/scripts/run_singlematch.py b/scripts/singlematch.py similarity index 100% rename from scripts/run_singlematch.py rename to scripts/singlematch.py diff --git a/scripts/run_split_halos.py b/scripts/split_halos.py similarity index 100% rename from scripts/run_split_halos.py rename to scripts/split_halos.py