From 4ee0a8d028a09722fca62d3db00d54a1ab0e810c Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Thu, 16 Feb 2023 12:57:24 +0100 Subject: [PATCH 01/14] Inserted a parameter so the size of outer components is not always increasing --- radiosim/jet.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index 4b99412..9b6af37 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -53,6 +53,7 @@ def create_jet(grid, conf): beta[1] = np.random.uniform(0, 1) y_rotation = np.random.uniform(0, np.pi) z_rotation = np.random.uniform(0, np.pi / 2) + extension = np.random.uniform(-0.5, 0.5) for i in range(comps): # amplitude decreases for more distant components, empirical @@ -89,8 +90,8 @@ def create_jet(grid, conf): img_size / comps * r_factor - * np.sqrt(i + 1) - / np.random.uniform(3, 8, size=2) + * (i + 1) ** extension + / np.random.uniform(5, 8, size=2) ) )[::-1] @@ -101,8 +102,8 @@ def create_jet(grid, conf): # print('Velocity of the jet:', beta) boost_app, boost_rec = relativistic_boosting(z_rotation, beta) - center_shift_x = np.random.uniform(-img_size / 20, img_size / 20) - center_shift_y = np.random.uniform(-img_size / 20, img_size / 20) + center_shift_x = np.random.uniform(-img_size / 50, img_size / 50) # will increase when zooming in + center_shift_y = np.random.uniform(-img_size / 50, img_size / 50) # will increase when zooming in if conf["scaling"] == "mojave": amp *= get_start_amp("mojave") From 51e53f607d8241026538514e2068caa57b8a4b44 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Wed, 22 Feb 2023 15:46:40 +0100 Subject: [PATCH 02/14] Change new parameters --- radiosim/jet.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index 9b6af37..0ddaf1e 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -53,7 +53,7 @@ def create_jet(grid, conf): beta[1] = np.random.uniform(0, 1) y_rotation = np.random.uniform(0, np.pi) z_rotation = np.random.uniform(0, np.pi / 2) - extension = np.random.uniform(-0.5, 0.5) + extension = np.random.uniform(0, 0.5) for i in range(comps): # amplitude decreases for more distant components, empirical @@ -91,7 +91,7 @@ def create_jet(grid, conf): / comps * r_factor * (i + 1) ** extension - / np.random.uniform(5, 8, size=2) + / np.random.uniform(3, 6, size=2) ) )[::-1] From b62bf2bed33c87f3a8d9f4be9acdf0e9a0065a0f Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Wed, 1 Mar 2023 14:21:24 +0100 Subject: [PATCH 03/14] Add posibility to save as fits file --- radiosim/jet.py | 2 +- radiosim/simulations.py | 5 ++-- radiosim/utils.py | 49 ++++++++++++++++++++++++++++++++++---- rc/default_simulation.toml | 1 + tests/simulate.toml | 1 + 5 files changed, 50 insertions(+), 8 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index 0ddaf1e..a8fb7b7 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -26,7 +26,7 @@ def create_jet(grid, conf): components. A jet with counter jet has c*2-1 components, since the center appears only once. Adding one channel for the backgound gives c*2 channels. source_lists: ndarray - array which stores all (seven) properties of each component, shape: [n, c*2-1, 7] + array which stores all (seven) properties of each component, shape: [n, c*2-1, 8] """ num_comps = conf["num_jet_components"] diff --git a/radiosim/simulations.py b/radiosim/simulations.py index 6a662fb..07e11c5 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -3,7 +3,6 @@ from radiosim.utils import ( create_grid, add_noise, - adjust_outpath, save_sky_distribution_bundle, ) from radiosim.jet import create_jet @@ -35,5 +34,5 @@ def create_sky_distribution(conf, opt: str): for img in sky_bundle: img -= img.min() img /= img.max() - path = adjust_outpath(conf["outpath"], "/samp_" + opt) - save_sky_distribution_bundle(path, sky_bundle, target_bundle) + + save_sky_distribution_bundle(conf, opt, sky_bundle, target_bundle) diff --git a/radiosim/utils.py b/radiosim/utils.py index f5ef011..7eed38d 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -9,6 +9,7 @@ from pathlib import Path from scipy import signal from astropy.convolution import Gaussian2DKernel +from astropy.io import fits def create_grid(pixel, bundle_size): @@ -165,7 +166,10 @@ def check_outpath(outpath, quiet=False): path = Path(outpath) exists = path.exists() if exists is True: - source = {p for p in path.rglob("*samp*.h5") if p.is_file()} + formats = ['.h5', '.fits'] + source = [] + for form in formats: + source.extend(p for p in path.rglob("*samp*" + form) if p.is_file()) if source: click.echo("Found existing source simulations!") if quiet: @@ -207,6 +211,7 @@ def read_config(config): """ sim_conf = {} sim_conf["mode"] = config["general"]["mode"] + sim_conf["fits"] = config["general"]["fits"] sim_conf["outpath"] = config["paths"]["outpath"] sim_conf["training_type"] = config["jet"]["training_type"] sim_conf["num_jet_components"] = config["jet"]["num_jet_components"] @@ -303,14 +308,16 @@ def adjust_outpath(path, option, form="h5"): return out -def save_sky_distribution_bundle(path, x, y, name_x="x", name_y="y"): +def save_sky_distribution_bundle(conf, opt, x, y, name_x="x", name_y="y"): """ Write images created in analysis to h5 file. Parameters ---------- - path: str - path to save file + conf: dict + configurations with all parameters + opt: str + dataset option (train, valid, test) x: ndarray image of the full jet, sum over all components y: ndarray @@ -320,11 +327,45 @@ def save_sky_distribution_bundle(path, x, y, name_x="x", name_y="y"): name_y: str name of the y-data """ + path = adjust_outpath(conf["outpath"], "/samp_" + opt) + with h5py.File(path, "w") as hf: hf.create_dataset(name_x, data=x) hf.create_dataset(name_y, data=y) hf.close() + if conf["fits"] and opt == "test": + for i in range(len(x)): + path_fits = adjust_outpath(conf["outpath"], "/samp_" + opt, form="fits") + + hdu_x = fits.PrimaryHDU(x[i, 0]) + hdul = fits.HDUList([hdu_x]) + + if conf["training_type"] == "list": + c1 = fits.Column(name='amplitude', array=y[i, :, 0], format='E') + c2 = fits.Column(name='x', array=y[i, :, 1], format='E') + c3 = fits.Column(name='y', array=y[i, :, 2], format='E') + c4 = fits.Column(name='sx', array=y[i, :, 3], format='E') + c5 = fits.Column(name='sy', array=y[i, :, 4], format='K') + c6 = fits.Column(name='rotation', array=y[i, :, 5], format='E') + c7 = fits.Column(name='z_rotation', array=y[i, :, 6], format='E') + c8 = fits.Column(name='beta', array=y[i, :, 7], format='E') + hdu_y = fits.BinTableHDU.from_columns([c1, c2, c3, c4, c5, c6, c7, c8]) + if y.shape[-1] != len(hdu_y.columns): + print(f"Simulated columns: {y.shape[-1]} vs. Saved columns: {len(hdu_y.columns)}") + hdul.append(hdu_y) + + elif conf["training_type"] == "clean": + hdu_y = fits.ImageHDU(y[i, 0], name=name_y) + hdul.append(hdu_y) + + elif conf["training_type"] == "gauss": + for j in range(y.shape[1]): + hdu_y = fits.ImageHDU(y[i, j], name=name_y + str(j)) + hdul.append(hdu_y) + + hdul.writeto(path_fits, overwrite=True) + def cart2pol(x: float, y: float): """ diff --git a/rc/default_simulation.toml b/rc/default_simulation.toml index ad1b164..4a68145 100644 --- a/rc/default_simulation.toml +++ b/rc/default_simulation.toml @@ -5,6 +5,7 @@ title = "Simulation configuration" [general] quiet = true mode = "jet" # jet or survey +fits = true [paths] outpath = "./build/example_data/" diff --git a/tests/simulate.toml b/tests/simulate.toml index 2ac01bc..5480460 100644 --- a/tests/simulate.toml +++ b/tests/simulate.toml @@ -5,6 +5,7 @@ title = "Simulation configuration" [general] quiet = true mode = "survey" # jet or survey +fits = true [paths] outpath = "./tests/build/data/" From b0f22d86620c6f314c52c45e8de7ccd9a77485d9 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Mon, 13 Mar 2023 17:09:05 +0100 Subject: [PATCH 04/14] Adjust header to PyBDSF demands --- radiosim/utils.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/radiosim/utils.py b/radiosim/utils.py index 7eed38d..37288c0 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -166,7 +166,7 @@ def check_outpath(outpath, quiet=False): path = Path(outpath) exists = path.exists() if exists is True: - formats = ['.h5', '.fits'] + formats = [".h5", ".fits"] source = [] for form in formats: source.extend(p for p in path.rglob("*samp*" + form) if p.is_file()) @@ -335,21 +335,28 @@ def save_sky_distribution_bundle(conf, opt, x, y, name_x="x", name_y="y"): hf.close() if conf["fits"] and opt == "test": + # header attributes needed for pybdsf + hdr = fits.Header() + hdr["CTYPE1"] = ("RA---SIN", "Axis name") + hdr["CDELT1"] = (-2.7777778241006E-08, "Pixel increment") + hdr["CTYPE2"] = ("DEC--SIN", "Axis name") + hdr["CDELT2"] = (2.7777778241006E-08, "Pixel increment") + for i in range(len(x)): path_fits = adjust_outpath(conf["outpath"], "/samp_" + opt, form="fits") - hdu_x = fits.PrimaryHDU(x[i, 0]) + hdu_x = fits.PrimaryHDU(data=x[i, 0], header=hdr) hdul = fits.HDUList([hdu_x]) if conf["training_type"] == "list": - c1 = fits.Column(name='amplitude', array=y[i, :, 0], format='E') - c2 = fits.Column(name='x', array=y[i, :, 1], format='E') - c3 = fits.Column(name='y', array=y[i, :, 2], format='E') - c4 = fits.Column(name='sx', array=y[i, :, 3], format='E') - c5 = fits.Column(name='sy', array=y[i, :, 4], format='K') - c6 = fits.Column(name='rotation', array=y[i, :, 5], format='E') - c7 = fits.Column(name='z_rotation', array=y[i, :, 6], format='E') - c8 = fits.Column(name='beta', array=y[i, :, 7], format='E') + c1 = fits.Column(name="amplitude", array=y[i, :, 0], format="E") + c2 = fits.Column(name="x", array=y[i, :, 1], format="E") + c3 = fits.Column(name="y", array=y[i, :, 2], format="E") + c4 = fits.Column(name="sx", array=y[i, :, 3], format="E") + c5 = fits.Column(name="sy", array=y[i, :, 4], format="K") + c6 = fits.Column(name="rotation", array=y[i, :, 5], format="E") + c7 = fits.Column(name="z_rotation", array=y[i, :, 6], format="E") + c8 = fits.Column(name="beta", array=y[i, :, 7], format="E") hdu_y = fits.BinTableHDU.from_columns([c1, c2, c3, c4, c5, c6, c7, c8]) if y.shape[-1] != len(hdu_y.columns): print(f"Simulated columns: {y.shape[-1]} vs. Saved columns: {len(hdu_y.columns)}") From 13881c5b2530ea122a4cb080b56c1d59c55f0a48 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Sun, 9 Apr 2023 10:16:25 +0200 Subject: [PATCH 05/14] Revert "Adjust header to PyBDSF demands" This reverts commit b0f22d86620c6f314c52c45e8de7ccd9a77485d9. --- radiosim/utils.py | 27 ++++++++++----------------- 1 file changed, 10 insertions(+), 17 deletions(-) diff --git a/radiosim/utils.py b/radiosim/utils.py index 37288c0..7eed38d 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -166,7 +166,7 @@ def check_outpath(outpath, quiet=False): path = Path(outpath) exists = path.exists() if exists is True: - formats = [".h5", ".fits"] + formats = ['.h5', '.fits'] source = [] for form in formats: source.extend(p for p in path.rglob("*samp*" + form) if p.is_file()) @@ -335,28 +335,21 @@ def save_sky_distribution_bundle(conf, opt, x, y, name_x="x", name_y="y"): hf.close() if conf["fits"] and opt == "test": - # header attributes needed for pybdsf - hdr = fits.Header() - hdr["CTYPE1"] = ("RA---SIN", "Axis name") - hdr["CDELT1"] = (-2.7777778241006E-08, "Pixel increment") - hdr["CTYPE2"] = ("DEC--SIN", "Axis name") - hdr["CDELT2"] = (2.7777778241006E-08, "Pixel increment") - for i in range(len(x)): path_fits = adjust_outpath(conf["outpath"], "/samp_" + opt, form="fits") - hdu_x = fits.PrimaryHDU(data=x[i, 0], header=hdr) + hdu_x = fits.PrimaryHDU(x[i, 0]) hdul = fits.HDUList([hdu_x]) if conf["training_type"] == "list": - c1 = fits.Column(name="amplitude", array=y[i, :, 0], format="E") - c2 = fits.Column(name="x", array=y[i, :, 1], format="E") - c3 = fits.Column(name="y", array=y[i, :, 2], format="E") - c4 = fits.Column(name="sx", array=y[i, :, 3], format="E") - c5 = fits.Column(name="sy", array=y[i, :, 4], format="K") - c6 = fits.Column(name="rotation", array=y[i, :, 5], format="E") - c7 = fits.Column(name="z_rotation", array=y[i, :, 6], format="E") - c8 = fits.Column(name="beta", array=y[i, :, 7], format="E") + c1 = fits.Column(name='amplitude', array=y[i, :, 0], format='E') + c2 = fits.Column(name='x', array=y[i, :, 1], format='E') + c3 = fits.Column(name='y', array=y[i, :, 2], format='E') + c4 = fits.Column(name='sx', array=y[i, :, 3], format='E') + c5 = fits.Column(name='sy', array=y[i, :, 4], format='K') + c6 = fits.Column(name='rotation', array=y[i, :, 5], format='E') + c7 = fits.Column(name='z_rotation', array=y[i, :, 6], format='E') + c8 = fits.Column(name='beta', array=y[i, :, 7], format='E') hdu_y = fits.BinTableHDU.from_columns([c1, c2, c3, c4, c5, c6, c7, c8]) if y.shape[-1] != len(hdu_y.columns): print(f"Simulated columns: {y.shape[-1]} vs. Saved columns: {len(hdu_y.columns)}") From 79242edce95e51b79e3b883f48afc22e1672f49b Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Sun, 9 Apr 2023 10:16:57 +0200 Subject: [PATCH 06/14] Revert "Add posibility to save as fits file" This reverts commit b62bf2bed33c87f3a8d9f4be9acdf0e9a0065a0f. --- radiosim/jet.py | 2 +- radiosim/simulations.py | 5 ++-- radiosim/utils.py | 49 ++++---------------------------------- rc/default_simulation.toml | 1 - tests/simulate.toml | 1 - 5 files changed, 8 insertions(+), 50 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index a8fb7b7..0ddaf1e 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -26,7 +26,7 @@ def create_jet(grid, conf): components. A jet with counter jet has c*2-1 components, since the center appears only once. Adding one channel for the backgound gives c*2 channels. source_lists: ndarray - array which stores all (seven) properties of each component, shape: [n, c*2-1, 8] + array which stores all (seven) properties of each component, shape: [n, c*2-1, 7] """ num_comps = conf["num_jet_components"] diff --git a/radiosim/simulations.py b/radiosim/simulations.py index 07e11c5..6a662fb 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -3,6 +3,7 @@ from radiosim.utils import ( create_grid, add_noise, + adjust_outpath, save_sky_distribution_bundle, ) from radiosim.jet import create_jet @@ -34,5 +35,5 @@ def create_sky_distribution(conf, opt: str): for img in sky_bundle: img -= img.min() img /= img.max() - - save_sky_distribution_bundle(conf, opt, sky_bundle, target_bundle) + path = adjust_outpath(conf["outpath"], "/samp_" + opt) + save_sky_distribution_bundle(path, sky_bundle, target_bundle) diff --git a/radiosim/utils.py b/radiosim/utils.py index 7eed38d..f5ef011 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -9,7 +9,6 @@ from pathlib import Path from scipy import signal from astropy.convolution import Gaussian2DKernel -from astropy.io import fits def create_grid(pixel, bundle_size): @@ -166,10 +165,7 @@ def check_outpath(outpath, quiet=False): path = Path(outpath) exists = path.exists() if exists is True: - formats = ['.h5', '.fits'] - source = [] - for form in formats: - source.extend(p for p in path.rglob("*samp*" + form) if p.is_file()) + source = {p for p in path.rglob("*samp*.h5") if p.is_file()} if source: click.echo("Found existing source simulations!") if quiet: @@ -211,7 +207,6 @@ def read_config(config): """ sim_conf = {} sim_conf["mode"] = config["general"]["mode"] - sim_conf["fits"] = config["general"]["fits"] sim_conf["outpath"] = config["paths"]["outpath"] sim_conf["training_type"] = config["jet"]["training_type"] sim_conf["num_jet_components"] = config["jet"]["num_jet_components"] @@ -308,16 +303,14 @@ def adjust_outpath(path, option, form="h5"): return out -def save_sky_distribution_bundle(conf, opt, x, y, name_x="x", name_y="y"): +def save_sky_distribution_bundle(path, x, y, name_x="x", name_y="y"): """ Write images created in analysis to h5 file. Parameters ---------- - conf: dict - configurations with all parameters - opt: str - dataset option (train, valid, test) + path: str + path to save file x: ndarray image of the full jet, sum over all components y: ndarray @@ -327,45 +320,11 @@ def save_sky_distribution_bundle(conf, opt, x, y, name_x="x", name_y="y"): name_y: str name of the y-data """ - path = adjust_outpath(conf["outpath"], "/samp_" + opt) - with h5py.File(path, "w") as hf: hf.create_dataset(name_x, data=x) hf.create_dataset(name_y, data=y) hf.close() - if conf["fits"] and opt == "test": - for i in range(len(x)): - path_fits = adjust_outpath(conf["outpath"], "/samp_" + opt, form="fits") - - hdu_x = fits.PrimaryHDU(x[i, 0]) - hdul = fits.HDUList([hdu_x]) - - if conf["training_type"] == "list": - c1 = fits.Column(name='amplitude', array=y[i, :, 0], format='E') - c2 = fits.Column(name='x', array=y[i, :, 1], format='E') - c3 = fits.Column(name='y', array=y[i, :, 2], format='E') - c4 = fits.Column(name='sx', array=y[i, :, 3], format='E') - c5 = fits.Column(name='sy', array=y[i, :, 4], format='K') - c6 = fits.Column(name='rotation', array=y[i, :, 5], format='E') - c7 = fits.Column(name='z_rotation', array=y[i, :, 6], format='E') - c8 = fits.Column(name='beta', array=y[i, :, 7], format='E') - hdu_y = fits.BinTableHDU.from_columns([c1, c2, c3, c4, c5, c6, c7, c8]) - if y.shape[-1] != len(hdu_y.columns): - print(f"Simulated columns: {y.shape[-1]} vs. Saved columns: {len(hdu_y.columns)}") - hdul.append(hdu_y) - - elif conf["training_type"] == "clean": - hdu_y = fits.ImageHDU(y[i, 0], name=name_y) - hdul.append(hdu_y) - - elif conf["training_type"] == "gauss": - for j in range(y.shape[1]): - hdu_y = fits.ImageHDU(y[i, j], name=name_y + str(j)) - hdul.append(hdu_y) - - hdul.writeto(path_fits, overwrite=True) - def cart2pol(x: float, y: float): """ diff --git a/rc/default_simulation.toml b/rc/default_simulation.toml index 4a68145..ad1b164 100644 --- a/rc/default_simulation.toml +++ b/rc/default_simulation.toml @@ -5,7 +5,6 @@ title = "Simulation configuration" [general] quiet = true mode = "jet" # jet or survey -fits = true [paths] outpath = "./build/example_data/" diff --git a/tests/simulate.toml b/tests/simulate.toml index 5480460..2ac01bc 100644 --- a/tests/simulate.toml +++ b/tests/simulate.toml @@ -5,7 +5,6 @@ title = "Simulation configuration" [general] quiet = true mode = "survey" # jet or survey -fits = true [paths] outpath = "./tests/build/data/" From 2c5884ea78723ce4daab37ace1a5fedab5e6de8d Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Sun, 9 Apr 2023 17:39:27 +0200 Subject: [PATCH 07/14] add class for multiprocessing --- radiosim/simulations.py | 43 ++++++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/radiosim/simulations.py b/radiosim/simulations.py index 6a662fb..bb69a60 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -1,4 +1,5 @@ import click +import multiprocessing from tqdm import tqdm from radiosim.utils import ( create_grid, @@ -12,28 +13,48 @@ def simulate_sky_distributions(conf): for opt in ["train", "valid", "test"]: - create_sky_distribution( + csd = create_sky_distribution( conf=conf, opt=opt, ) + csd() -def create_sky_distribution(conf, opt: str): - for _ in tqdm(range(conf["bundles_" + opt])): - grid = create_grid(conf["img_size"], conf["bundle_size"]) - if conf["mode"] == "jet": - sky, target = create_jet(grid, conf) - elif conf["mode"] == "survey": - sky, target = create_survey(grid, conf) +class create_sky_distribution: + def __init__(self, conf, opt): + self.conf = conf + self.opt = opt + + def __call__(self): + n_bundels = self.conf["bundles_" + self.opt] + n_cores = int(multiprocessing.cpu_count() * 0.5) # use 50% of available cores + print("Number of cpu cores:", n_cores) + + if n_cores == 1: + for _ in tqdm(range(n_bundels)): + self.multiprocessing_func(0) + else: + print() + with multiprocessing.Pool(n_cores) as p: + # r = list(tqdm(p.imap(self.multiprocessing_func, range(n_bundels)), total=n_bundels)) # sometimes leads to error in tqdm: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') + for _ in p.imap(self.multiprocessing_func, range(n_bundels)): + continue + + def multiprocessing_func(self, _): + grid = create_grid(self.conf["img_size"], self.conf["bundle_size"]) + if self.conf["mode"] == "jet": + sky, target = create_jet(grid, self.conf) + elif self.conf["mode"] == "survey": + sky, target = create_survey(grid, self.conf) else: click.echo("Given mode not found. Choose 'survey' or 'jet' in config file") sky_bundle = sky.copy() target_bundle = target.copy() - if conf["noise"] and conf["noise_level"] > 0: - sky_bundle = add_noise(sky_bundle, conf["noise_level"]) + if self.conf["noise"] and self.conf["noise_level"] > 0: + sky_bundle = add_noise(sky_bundle, self.conf["noise_level"]) for img in sky_bundle: img -= img.min() img /= img.max() - path = adjust_outpath(conf["outpath"], "/samp_" + opt) + path = adjust_outpath(self.conf["outpath"], "/samp_" + self.opt) save_sky_distribution_bundle(path, sky_bundle, target_bundle) From 49aa1ad71bfb4e42c22156e6cb33f092097c7e01 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Tue, 11 Apr 2023 08:38:25 +0200 Subject: [PATCH 08/14] added to config file --- radiosim/scripts/start_simulation.py | 4 ++-- radiosim/simulations.py | 19 ++++++++++++------- radiosim/utils.py | 2 ++ rc/default_simulation.toml | 3 ++- tests/simulate.toml | 3 ++- 5 files changed, 20 insertions(+), 11 deletions(-) diff --git a/radiosim/scripts/start_simulation.py b/radiosim/scripts/start_simulation.py index 488fa9e..045b578 100644 --- a/radiosim/scripts/start_simulation.py +++ b/radiosim/scripts/start_simulation.py @@ -22,10 +22,10 @@ def main(configuration_path, mode): conf = read_config(config) print(conf, "\n") - outpath = config["paths"]["outpath"] + outpath = conf["outpath"] sim_sources = check_outpath( outpath, - quiet=config["general"]["quiet"], + quiet=conf["quiet"], ) if mode == "simulate": diff --git a/radiosim/simulations.py b/radiosim/simulations.py index bb69a60..3dbefb4 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -28,19 +28,24 @@ def __init__(self, conf, opt): def __call__(self): n_bundels = self.conf["bundles_" + self.opt] n_cores = int(multiprocessing.cpu_count() * 0.5) # use 50% of available cores - print("Number of cpu cores:", n_cores) - - if n_cores == 1: + if n_cores == 1 or not self.conf["multiprocessing"]: for _ in tqdm(range(n_bundels)): - self.multiprocessing_func(0) + self.sky_distribution(0) else: print() with multiprocessing.Pool(n_cores) as p: - # r = list(tqdm(p.imap(self.multiprocessing_func, range(n_bundels)), total=n_bundels)) # sometimes leads to error in tqdm: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') - for _ in p.imap(self.multiprocessing_func, range(n_bundels)): + # r = list(tqdm(p.imap(self.sky_distribution, range(n_bundels)), total=n_bundels)) # sometimes leads to error in tqdm: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') + for _ in p.imap(self.sky_distribution, range(n_bundels)): continue - def multiprocessing_func(self, _): + def sky_distribution(self, _): + """Create and save the sky distribution + + Parameters + ---------- + _: Any + dummy parameter needed for 'imap' method of multiprocessing + """ grid = create_grid(self.conf["img_size"], self.conf["bundle_size"]) if self.conf["mode"] == "jet": sky, target = create_jet(grid, self.conf) diff --git a/radiosim/utils.py b/radiosim/utils.py index f5ef011..d0a1d7f 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -206,7 +206,9 @@ def read_config(config): unpacked configurations """ sim_conf = {} + sim_conf["quiet"] = config["general"]["quiet"] sim_conf["mode"] = config["general"]["mode"] + sim_conf["multiprocessing"] = config["general"]["multiprocessing"] sim_conf["outpath"] = config["paths"]["outpath"] sim_conf["training_type"] = config["jet"]["training_type"] sim_conf["num_jet_components"] = config["jet"]["num_jet_components"] diff --git a/rc/default_simulation.toml b/rc/default_simulation.toml index ad1b164..4129d4d 100644 --- a/rc/default_simulation.toml +++ b/rc/default_simulation.toml @@ -4,7 +4,8 @@ title = "Simulation configuration" [general] quiet = true -mode = "jet" # jet or survey +mode = "jet" # jet or survey +multiprocessing = true # no progress bar (to be implemented) [paths] outpath = "./build/example_data/" diff --git a/tests/simulate.toml b/tests/simulate.toml index 2ac01bc..ec5ab30 100644 --- a/tests/simulate.toml +++ b/tests/simulate.toml @@ -4,7 +4,8 @@ title = "Simulation configuration" [general] quiet = true -mode = "survey" # jet or survey +mode = "survey" # jet or survey +multiprocessing = true # no progress bar (to be implemented) [paths] outpath = "./tests/build/data/" From 86dbc26ba03c490f35f78739e760173a5bd0be58 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Mon, 17 Apr 2023 21:26:30 +0200 Subject: [PATCH 09/14] Adjust parameters, fix y-rotation (only 180 degree) --- radiosim/jet.py | 64 ++++++++++++++++++++--------------------- radiosim/simulations.py | 8 ++++-- rc/simulation.toml | 31 ++++++++++++++++++++ 3 files changed, 68 insertions(+), 35 deletions(-) create mode 100644 rc/simulation.toml diff --git a/radiosim/jet.py b/radiosim/jet.py index 0ddaf1e..0623163 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -51,7 +51,7 @@ def create_jet(grid, conf): # velocity in units of c_0, initialise velocity of first component beta = np.zeros(num_comps[1]) beta[1] = np.random.uniform(0, 1) - y_rotation = np.random.uniform(0, np.pi) + y_rotation = np.random.uniform(0, 2 * np.pi) z_rotation = np.random.uniform(0, np.pi / 2) extension = np.random.uniform(0, 0.5) @@ -66,7 +66,7 @@ def create_jet(grid, conf): beta[i] = beta[1] * np.exp(-np.sqrt(i - 1) * np.random.normal(0.5, 0.1)) # curving the jet, empirical - y_rotation += np.random.normal(0, np.pi / 24) + y_rotation += np.random.normal(0, np.pi / 36) # distance between components, r_factor to fill the corners jet_angle_cos = np.abs(np.cos(y_rotation)) @@ -91,7 +91,7 @@ def create_jet(grid, conf): / comps * r_factor * (i + 1) ** extension - / np.random.uniform(3, 6, size=2) + / np.random.uniform(3, 9, size=2) ) )[::-1] @@ -102,8 +102,8 @@ def create_jet(grid, conf): # print('Velocity of the jet:', beta) boost_app, boost_rec = relativistic_boosting(z_rotation, beta) - center_shift_x = np.random.uniform(-img_size / 50, img_size / 50) # will increase when zooming in - center_shift_y = np.random.uniform(-img_size / 50, img_size / 50) # will increase when zooming in + center_shift_x = np.random.uniform(-img_size / 20, img_size / 20) # will increase when zooming in + center_shift_y = np.random.uniform(-img_size / 20, img_size / 20) # will increase when zooming in if conf["scaling"] == "mojave": amp *= get_start_amp("mojave") @@ -133,33 +133,33 @@ def create_jet(grid, conf): jet_comp = np.array(component_from_list(img_size, amp, x, y, sx, sy, rotation)) jet_img = np.sum(jet_comp, axis=0) - # get values at the edge of the image - edge_list = [ - jet_img[0, :-1], - jet_img[:-1, -1], - jet_img[-1, ::-1], - jet_img[-2:0:-1, 0], - ] - edges = np.concatenate(edge_list) - edge_threshold = 0.01 - if edges.max() < edge_threshold: - # zoom on source to equalize size differences from z-rotation - jet_img, jet_comp, zoom_factor = zoom_on_source( - jet_img, jet_comp, max_amp=edge_threshold - ) - x = img_size / 2 + (x - img_size / 2) * zoom_factor - y = img_size / 2 + (y - img_size / 2) * zoom_factor - sx *= zoom_factor - sy *= zoom_factor - - # random zoom out for more variance - zoom_out_factor = np.random.uniform(1 / 2, 1) # 1/8: pad eg. 128 -> 1024 - pad_value = (1 / zoom_out_factor - 1) * img_size / 2 - jet_img, jet_comp = zoom_out(jet_img, jet_comp, pad_value=pad_value) - x = img_size / 2 + (x - img_size / 2) * zoom_out_factor - y = img_size / 2 + (y - img_size / 2) * zoom_out_factor - sx *= zoom_out_factor - sy *= zoom_out_factor + # # get values at the edge of the image + # edge_list = [ + # jet_img[0, :-1], + # jet_img[:-1, -1], + # jet_img[-1, ::-1], + # jet_img[-2:0:-1, 0], + # ] + # edges = np.concatenate(edge_list) + # edge_threshold = 0.01 + # if edges.max() < edge_threshold: + # # zoom on source to equalize size differences from z-rotation + # jet_img, jet_comp, zoom_factor = zoom_on_source( + # jet_img, jet_comp, max_amp=edge_threshold + # ) + # x = img_size / 2 + (x - img_size / 2) * zoom_factor + # y = img_size / 2 + (y - img_size / 2) * zoom_factor + # sx *= zoom_factor + # sy *= zoom_factor + + # # random zoom out for more variance + # zoom_out_factor = np.random.uniform(1 / 2, 1) # 1/8: pad eg. 128 -> 1024 + # pad_value = (1 / zoom_out_factor - 1) * img_size / 2 + # jet_img, jet_comp = zoom_out(jet_img, jet_comp, pad_value=pad_value) + # x = img_size / 2 + (x - img_size / 2) * zoom_out_factor + # y = img_size / 2 + (y - img_size / 2) * zoom_out_factor + # sx *= zoom_out_factor + # sy *= zoom_out_factor # normalisation if conf["scaling"] == "normalize": diff --git a/radiosim/simulations.py b/radiosim/simulations.py index 3dbefb4..6208958 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -34,9 +34,10 @@ def __call__(self): else: print() with multiprocessing.Pool(n_cores) as p: - # r = list(tqdm(p.imap(self.sky_distribution, range(n_bundels)), total=n_bundels)) # sometimes leads to error in tqdm: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') - for _ in p.imap(self.sky_distribution, range(n_bundels)): - continue + _ = list(tqdm(p.imap(self.sky_distribution, range(n_bundels)), total=n_bundels)) # sometimes leads to error in tqdm: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') + + # for _ in p.imap(self.sky_distribution, range(n_bundels)): + # continue def sky_distribution(self, _): """Create and save the sky distribution @@ -63,3 +64,4 @@ def sky_distribution(self, _): img /= img.max() path = adjust_outpath(self.conf["outpath"], "/samp_" + self.opt) save_sky_distribution_bundle(path, sky_bundle, target_bundle) + \ No newline at end of file diff --git a/rc/simulation.toml b/rc/simulation.toml new file mode 100644 index 0000000..d3a60e2 --- /dev/null +++ b/rc/simulation.toml @@ -0,0 +1,31 @@ +# This is a TOML document. + +title = "Simulation configuration" + +[general] +quiet = true +mode = "jet" # jet or survey +multiprocessing = true # no progress bar (to be implemented) + +[paths] +outpath = "/mnt/nvme1n1/simulations/230417_n10000_s256" + +[jet] +training_type = "list" # gauss or list or clean +num_jet_components = [3, 9] +scaling = "normalize" + +[survey] +num_sources = 20 +# classes are: jet, gaussian, pointlike +class_distribution = [2, 1, 2] # average share of each source in the images +scale_sources = true + +[image_options] +bundles_train = 100 +bundles_valid = 20 +bundles_test = 20 +bundle_size = 100 +img_size = 256 +noise = true +noise_level = 15 # maximum intensity of noise in percent From 822e953ab2ce9e67de89c0ac8aaddc057f5497f3 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Mon, 24 Apr 2023 10:06:45 +0200 Subject: [PATCH 10/14] Changed droprate and multiprocessing --- radiosim/jet.py | 4 ++-- radiosim/simulations.py | 17 +++++++++++++++-- rc/default_simulation.toml | 2 +- rc/simulation.toml | 10 +++++----- tests/simulate.toml | 2 +- 5 files changed, 24 insertions(+), 11 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index 0623163..9313604 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -95,7 +95,7 @@ def create_jet(grid, conf): ) )[::-1] - # rotation aligned with the jet angle, empirical + # rotation of the gauss component, empirical if i >= 1: rotation[i] = rotation[i - 1] + np.random.normal(0, np.pi / 18) @@ -110,7 +110,7 @@ def create_jet(grid, conf): # mirror the data for the counter jet # random drop of counter jet, because the relativistic boosting only does not create clear one-sided jets - if np.random.rand() < 0.3: + if np.random.rand() < 0.5: amp = np.concatenate((amp * boost_app, amp[1:] * boost_rec[1:])) x = np.concatenate((x + center_shift_x, img_size - x[1:] + center_shift_x)) y = np.concatenate((y + center_shift_y, img_size - y[1:] + center_shift_y)) diff --git a/radiosim/simulations.py b/radiosim/simulations.py index 6208958..0f0c31e 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -33,11 +33,24 @@ def __call__(self): self.sky_distribution(0) else: print() + # with multiprocessing.Pool(n_cores) as p: + # with tqdm(total=n_bundels) as pbar: + # for _ in p.imap(self.sky_distribution, range(n_bundels)): + # pbar.update() + # # p.close() + # # p.join() + with multiprocessing.Pool(n_cores) as p: - _ = list(tqdm(p.imap(self.sky_distribution, range(n_bundels)), total=n_bundels)) # sometimes leads to error in tqdm: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') + # can lead to error: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') + + # with tqdm + results = list(tqdm(p.imap(self.sky_distribution, range(n_bundels)), total=n_bundels)) + # for result in results: + # pass + # without tqdm # for _ in p.imap(self.sky_distribution, range(n_bundels)): - # continue + # pass def sky_distribution(self, _): """Create and save the sky distribution diff --git a/rc/default_simulation.toml b/rc/default_simulation.toml index 4129d4d..9c87fa7 100644 --- a/rc/default_simulation.toml +++ b/rc/default_simulation.toml @@ -5,7 +5,7 @@ title = "Simulation configuration" [general] quiet = true mode = "jet" # jet or survey -multiprocessing = true # no progress bar (to be implemented) +multiprocessing = true [paths] outpath = "./build/example_data/" diff --git a/rc/simulation.toml b/rc/simulation.toml index d3a60e2..be463d2 100644 --- a/rc/simulation.toml +++ b/rc/simulation.toml @@ -5,10 +5,10 @@ title = "Simulation configuration" [general] quiet = true mode = "jet" # jet or survey -multiprocessing = true # no progress bar (to be implemented) +multiprocessing = true [paths] -outpath = "/mnt/nvme1n1/simulations/230417_n10000_s256" +outpath = "/mnt/nvme1n1/simulations/230423_n50000_s256" [jet] training_type = "list" # gauss or list or clean @@ -22,9 +22,9 @@ class_distribution = [2, 1, 2] # average share of each source in the images scale_sources = true [image_options] -bundles_train = 100 -bundles_valid = 20 -bundles_test = 20 +bundles_train = 500 +bundles_valid = 50 +bundles_test = 50 bundle_size = 100 img_size = 256 noise = true diff --git a/tests/simulate.toml b/tests/simulate.toml index ec5ab30..bf41d33 100644 --- a/tests/simulate.toml +++ b/tests/simulate.toml @@ -5,7 +5,7 @@ title = "Simulation configuration" [general] quiet = true mode = "survey" # jet or survey -multiprocessing = true # no progress bar (to be implemented) +multiprocessing = true [paths] outpath = "./tests/build/data/" From 381a25d70487e76a0b300ece501ddabe47f9bd21 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Mon, 1 May 2023 10:39:21 +0200 Subject: [PATCH 11/14] Fix multiprocessing BlockingIOError --- radiosim/jet.py | 10 +++++--- radiosim/simulations.py | 46 +++++++++++++----------------------- radiosim/utils.py | 52 ++++++++++++++--------------------------- rc/simulation.toml | 6 ++--- 4 files changed, 44 insertions(+), 70 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index 9313604..bb9377e 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -57,7 +57,7 @@ def create_jet(grid, conf): for i in range(comps): # amplitude decreases for more distant components, empirical - amp[i] = np.exp(-np.sqrt(i) * np.random.normal(1.3, 0.2)) + amp[i] = np.exp(-np.sqrt(i) * np.random.normal(1.1, 0.2)) if i >= 1 and np.random.rand() < 0.1: # drop some components amp[i] = 0 @@ -102,8 +102,12 @@ def create_jet(grid, conf): # print('Velocity of the jet:', beta) boost_app, boost_rec = relativistic_boosting(z_rotation, beta) - center_shift_x = np.random.uniform(-img_size / 20, img_size / 20) # will increase when zooming in - center_shift_y = np.random.uniform(-img_size / 20, img_size / 20) # will increase when zooming in + center_shift_x = np.random.uniform( + -img_size / 20, img_size / 20 + ) # will increase when zooming in + center_shift_y = np.random.uniform( + -img_size / 20, img_size / 20 + ) # will increase when zooming in if conf["scaling"] == "mojave": amp *= get_start_amp("mojave") diff --git a/radiosim/simulations.py b/radiosim/simulations.py index 0f0c31e..20db4d5 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -1,14 +1,14 @@ import click import multiprocessing from tqdm import tqdm + +from radiosim.jet import create_jet +from radiosim.survey import create_survey from radiosim.utils import ( - create_grid, add_noise, - adjust_outpath, + create_grid, save_sky_distribution_bundle, ) -from radiosim.jet import create_jet -from radiosim.survey import create_survey def simulate_sky_distributions(conf): @@ -29,36 +29,23 @@ def __call__(self): n_bundels = self.conf["bundles_" + self.opt] n_cores = int(multiprocessing.cpu_count() * 0.5) # use 50% of available cores if n_cores == 1 or not self.conf["multiprocessing"]: - for _ in tqdm(range(n_bundels)): - self.sky_distribution(0) + for i in tqdm(range(n_bundels)): + self.sky_distribution(i) else: - print() - # with multiprocessing.Pool(n_cores) as p: - # with tqdm(total=n_bundels) as pbar: - # for _ in p.imap(self.sky_distribution, range(n_bundels)): - # pbar.update() - # # p.close() - # # p.join() - with multiprocessing.Pool(n_cores) as p: - # can lead to error: BlockingIOError: [Errno 11] Unable to create file (unable to lock file, errno = 11, error message = 'Resource temporarily unavailable') + _ = list( + tqdm( + p.imap(self.sky_distribution, range(n_bundels)), total=n_bundels + ) + ) - # with tqdm - results = list(tqdm(p.imap(self.sky_distribution, range(n_bundels)), total=n_bundels)) - # for result in results: - # pass - - # without tqdm - # for _ in p.imap(self.sky_distribution, range(n_bundels)): - # pass - - def sky_distribution(self, _): + def sky_distribution(self, i: int): """Create and save the sky distribution - + Parameters ---------- - _: Any - dummy parameter needed for 'imap' method of multiprocessing + i: int + n-th sky distribution to be saved """ grid = create_grid(self.conf["img_size"], self.conf["bundle_size"]) if self.conf["mode"] == "jet": @@ -75,6 +62,5 @@ def sky_distribution(self, _): for img in sky_bundle: img -= img.min() img /= img.max() - path = adjust_outpath(self.conf["outpath"], "/samp_" + self.opt) + path = self.conf["outpath"] + "/samp_" + self.opt + "_" + str(i) + ".h5" save_sky_distribution_bundle(path, sky_bundle, target_bundle) - \ No newline at end of file diff --git a/radiosim/utils.py b/radiosim/utils.py index d0a1d7f..dcf2254 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -57,7 +57,7 @@ def relativistic_boosting(theta, beta): boost_rec: float boosting factor for the receding jet """ - gamma = 1 / np.sqrt(1 - beta ** 2) # Lorentz factor + gamma = 1 / np.sqrt(1 - beta**2) # Lorentz factor mu = np.cos(theta) boost_app = 1 / (gamma * (1 - beta * mu)) @@ -99,14 +99,18 @@ def zoom_on_source(img, comp=None, max_amp=0.01): zoom_factor = size / (size - 2 * idx) # crop the source - cropped_img = img[idx:size-idx, idx:size-idx] - zoomed_img = cv2.resize(cropped_img, dsize=(size, size), interpolation=cv2.INTER_LINEAR) + cropped_img = img[idx : size - idx, idx : size - idx] + zoomed_img = cv2.resize( + cropped_img, dsize=(size, size), interpolation=cv2.INTER_LINEAR + ) if comp is not None: - cropped_comp = comp[:, idx:size-idx, idx:size-idx] + cropped_comp = comp[:, idx : size - idx, idx : size - idx] zoomed_comp = np.empty_like(comp) for i, component in enumerate(cropped_comp): - zoomed_comp[i] = cv2.resize(component, dsize=(size, size), interpolation=cv2.INTER_LINEAR) + zoomed_comp[i] = cv2.resize( + component, dsize=(size, size), interpolation=cv2.INTER_LINEAR + ) return zoomed_img, zoomed_comp, zoom_factor return zoomed_img, zoom_factor @@ -135,11 +139,17 @@ def zoom_out(img, comp=None, pad_value=0): if not isinstance(pad_value, int): pad_value = np.int64(pad_value) size = img.shape[0] - img = cv2.resize(np.pad(img, pad_value), dsize=(size, size), interpolation=cv2.INTER_LINEAR) + img = cv2.resize( + np.pad(img, pad_value), dsize=(size, size), interpolation=cv2.INTER_LINEAR + ) if comp is not None: for component in comp: - component = cv2.resize(np.pad(component, pad_value), dsize=(size, size), interpolation=cv2.INTER_LINEAR) + component = cv2.resize( + np.pad(component, pad_value), + dsize=(size, size), + interpolation=cv2.INTER_LINEAR, + ) return img, comp return img @@ -279,32 +289,6 @@ def call_noise(kernels, strengths): return image_noised -def adjust_outpath(path, option, form="h5"): - """ - Add number to out path when filename already exists. - - Parameters - ---------- - path: str - path to save directory - option: str - additional keyword to add to path - form: str - file extension - - Returns - ------- - out: str - adjusted path - """ - counter = 0 - filename = str(path) + (option + "_{}." + form) - while os.path.isfile(filename.format(counter)): - counter += 1 - out = filename.format(counter) - return out - - def save_sky_distribution_bundle(path, x, y, name_x="x", name_y="y"): """ Write images created in analysis to h5 file. @@ -346,7 +330,7 @@ def cart2pol(x: float, y: float): phi: float angle in radian """ - r = np.sqrt(x ** 2 + y ** 2) + r = np.sqrt(x**2 + y**2) phi = np.arctan2(y, x) return (r, phi) diff --git a/rc/simulation.toml b/rc/simulation.toml index be463d2..fecadb9 100644 --- a/rc/simulation.toml +++ b/rc/simulation.toml @@ -8,11 +8,11 @@ mode = "jet" # jet or survey multiprocessing = true [paths] -outpath = "/mnt/nvme1n1/simulations/230423_n50000_s256" +outpath = "/mnt/nvme1n1/simulations/230502_n50000_s128" [jet] training_type = "list" # gauss or list or clean -num_jet_components = [3, 9] +num_jet_components = [2, 9] scaling = "normalize" [survey] @@ -26,6 +26,6 @@ bundles_train = 500 bundles_valid = 50 bundles_test = 50 bundle_size = 100 -img_size = 256 +img_size = 128 noise = true noise_level = 15 # maximum intensity of noise in percent From bbf1d9e6dd4647c5a8254bc6bac0eb52b8faa0c1 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Mon, 1 May 2023 11:06:19 +0200 Subject: [PATCH 12/14] Renaming and deletion for better readability --- radiosim/jet.py | 34 ++---------------- radiosim/utils.py | 90 ----------------------------------------------- 2 files changed, 3 insertions(+), 121 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index bb9377e..ba4a07b 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -1,5 +1,5 @@ import numpy as np -from radiosim.utils import relativistic_boosting, pol2cart, zoom_on_source, zoom_out +from radiosim.utils import relativistic_boosting, pol2cart from radiosim.gauss import twodgaussian from radiosim.flux_scaling import get_start_amp @@ -53,7 +53,7 @@ def create_jet(grid, conf): beta[1] = np.random.uniform(0, 1) y_rotation = np.random.uniform(0, 2 * np.pi) z_rotation = np.random.uniform(0, np.pi / 2) - extension = np.random.uniform(0, 0.5) + expansion = np.random.uniform(0, 0.5) for i in range(comps): # amplitude decreases for more distant components, empirical @@ -90,7 +90,7 @@ def create_jet(grid, conf): img_size / comps * r_factor - * (i + 1) ** extension + * (i + 1) ** expansion / np.random.uniform(3, 9, size=2) ) )[::-1] @@ -137,34 +137,6 @@ def create_jet(grid, conf): jet_comp = np.array(component_from_list(img_size, amp, x, y, sx, sy, rotation)) jet_img = np.sum(jet_comp, axis=0) - # # get values at the edge of the image - # edge_list = [ - # jet_img[0, :-1], - # jet_img[:-1, -1], - # jet_img[-1, ::-1], - # jet_img[-2:0:-1, 0], - # ] - # edges = np.concatenate(edge_list) - # edge_threshold = 0.01 - # if edges.max() < edge_threshold: - # # zoom on source to equalize size differences from z-rotation - # jet_img, jet_comp, zoom_factor = zoom_on_source( - # jet_img, jet_comp, max_amp=edge_threshold - # ) - # x = img_size / 2 + (x - img_size / 2) * zoom_factor - # y = img_size / 2 + (y - img_size / 2) * zoom_factor - # sx *= zoom_factor - # sy *= zoom_factor - - # # random zoom out for more variance - # zoom_out_factor = np.random.uniform(1 / 2, 1) # 1/8: pad eg. 128 -> 1024 - # pad_value = (1 / zoom_out_factor - 1) * img_size / 2 - # jet_img, jet_comp = zoom_out(jet_img, jet_comp, pad_value=pad_value) - # x = img_size / 2 + (x - img_size / 2) * zoom_out_factor - # y = img_size / 2 + (y - img_size / 2) * zoom_out_factor - # sx *= zoom_out_factor - # sy *= zoom_out_factor - # normalisation if conf["scaling"] == "normalize": jet_max = jet_img.max() diff --git a/radiosim/utils.py b/radiosim/utils.py index dcf2254..ef6a431 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -65,96 +65,6 @@ def relativistic_boosting(theta, beta): return boost_app, boost_rec -def zoom_on_source(img, comp=None, max_amp=0.01): - """ - Zoom on source to cut out irrelevant area. Shape will stay equal. - - Parameters - ---------- - img: 2D array - Image of the sky used to zoom on - comp: 3D array (n, (img)) - Images of the components, same zooming as on img - max_amp: float - Maximal amplitude which will be at the edge of the image - - Returns - ------- - zoomed_img: ndarray - Image after zooming - zoom_factor: float - Zooming factor - """ - # find farest outside column or row with amplitude > max_amp - mask = img > max_amp - mask_flip = np.flip(mask) - - idx_left = np.argmax(np.sum(mask, axis=0) > 0) - idx_right = np.argmax(np.sum(mask_flip, axis=0) > 0) - idx_bottom = np.argmax(np.sum(mask, axis=1) > 0) - idx_top = np.argmax(np.sum(mask_flip, axis=1) > 0) - # print(idx_left, idx_right, idx_bottom, idx_top) - idx = np.min([idx_left, idx_right, idx_bottom, idx_top]) - size = img.shape[0] - zoom_factor = size / (size - 2 * idx) - - # crop the source - cropped_img = img[idx : size - idx, idx : size - idx] - zoomed_img = cv2.resize( - cropped_img, dsize=(size, size), interpolation=cv2.INTER_LINEAR - ) - - if comp is not None: - cropped_comp = comp[:, idx : size - idx, idx : size - idx] - zoomed_comp = np.empty_like(comp) - for i, component in enumerate(cropped_comp): - zoomed_comp[i] = cv2.resize( - component, dsize=(size, size), interpolation=cv2.INTER_LINEAR - ) - return zoomed_img, zoomed_comp, zoom_factor - - return zoomed_img, zoom_factor - - -def zoom_out(img, comp=None, pad_value=0): - """ - Zoom out of an image. Padding edges with zeros. - - Parameters - ---------- - img: 2D array - Image of the sky - comp: 3D array (n, (img)) - Images of the components - pad_value: int - Number of pixels to pad around the source - - Returns - ------- - zoomed_img: ndarray - Image after zooming - zoomed_comp: ndarray - Componets after zooming - """ - if not isinstance(pad_value, int): - pad_value = np.int64(pad_value) - size = img.shape[0] - img = cv2.resize( - np.pad(img, pad_value), dsize=(size, size), interpolation=cv2.INTER_LINEAR - ) - - if comp is not None: - for component in comp: - component = cv2.resize( - np.pad(component, pad_value), - dsize=(size, size), - interpolation=cv2.INTER_LINEAR, - ) - return img, comp - - return img - - def check_outpath(outpath, quiet=False): """ Check if outpath exists. Check for existing source_bundle files. From d8f89755ffb4aa6c4cadea866962110849d9c333 Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Mon, 1 May 2023 11:11:18 +0200 Subject: [PATCH 13/14] Remove own simulation toml --- rc/simulation.toml | 31 ------------------------------- 1 file changed, 31 deletions(-) delete mode 100644 rc/simulation.toml diff --git a/rc/simulation.toml b/rc/simulation.toml deleted file mode 100644 index fecadb9..0000000 --- a/rc/simulation.toml +++ /dev/null @@ -1,31 +0,0 @@ -# This is a TOML document. - -title = "Simulation configuration" - -[general] -quiet = true -mode = "jet" # jet or survey -multiprocessing = true - -[paths] -outpath = "/mnt/nvme1n1/simulations/230502_n50000_s128" - -[jet] -training_type = "list" # gauss or list or clean -num_jet_components = [2, 9] -scaling = "normalize" - -[survey] -num_sources = 20 -# classes are: jet, gaussian, pointlike -class_distribution = [2, 1, 2] # average share of each source in the images -scale_sources = true - -[image_options] -bundles_train = 500 -bundles_valid = 50 -bundles_test = 50 -bundle_size = 100 -img_size = 128 -noise = true -noise_level = 15 # maximum intensity of noise in percent From cb3a0e14e9f5b3e63b63d90894da20ec0292b5ff Mon Sep 17 00:00:00 2001 From: Arne Poggenpohl Date: Thu, 4 May 2023 11:01:48 +0200 Subject: [PATCH 14/14] Number of components as limit for position, not as fixed value --- radiosim/jet.py | 9 +++++---- setup.py | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/radiosim/jet.py b/radiosim/jet.py index ba4a07b..3f3689f 100644 --- a/radiosim/jet.py +++ b/radiosim/jet.py @@ -39,6 +39,7 @@ def create_jet(grid, conf): targets = [] for _ in grid: comps = np.random.randint(num_comps[0], num_comps[1] + 1) + comps4sim = np.random.uniform(comps, num_comps[1]) amp = np.zeros(num_comps[1]) x = np.zeros(num_comps[1]) @@ -79,7 +80,7 @@ def create_jet(grid, conf): r_factor = np.sqrt(2) # *0.7 so the last component is not on the edge - r = i / (comps - 1) * img_size / 2 * r_factor * np.sin(z_rotation) * 0.7 + r = i / (comps4sim - 1) * img_size / 2 * r_factor * np.sin(z_rotation) * 0.7 # get the cartesian coordinates x[i], y[i] = np.array(pol2cart(r, y_rotation)) + center @@ -88,7 +89,7 @@ def create_jet(grid, conf): sx[i], sy[i] = np.sort( ( img_size - / comps + / comps4sim * r_factor * (i + 1) ** expansion / np.random.uniform(3, 9, size=2) @@ -104,10 +105,10 @@ def create_jet(grid, conf): center_shift_x = np.random.uniform( -img_size / 20, img_size / 20 - ) # will increase when zooming in + ) center_shift_y = np.random.uniform( -img_size / 20, img_size / 20 - ) # will increase when zooming in + ) if conf["scaling"] == "mojave": amp *= get_start_amp("mojave") diff --git a/setup.py b/setup.py index 499d745..d5e43c3 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ setup( name="radiosim", - version="0.0.2", + version="0.0.3", description="Simulation of radio skies to create astrophysical data sets", url="https://github.com/Kevin2/radionets", author="Kevin Schmidt, Felix Geyer, Paul-Simon Blomenkamp, Stefan Fröse",