diff --git a/README.md b/README.md index b68303c..029487b 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,30 @@ # radiosim Simulation of radio skies to create astrophysical data sets. +This repository is part of the [`radionets-project`](https://github.com/radionets-project). + +## Installation + +This repository is built as a python package. We recommend creating a mamba environment to handle the dependencies of all packages. +You can create one by running the following command in this repository: +``` +$ mamba env create -f environment.yml +``` +You can start the environment with +``` +$ mamba activate radiosim +``` +after the installation. + +## Usage + +There are currently three supported simulation types: +1. `survey` full sky simulation +2. `jet` extended source +3. `mojave` MOJAVE like extended source + +In the `radiosim` environment you can start the simulation with +``` +$ radiosim path/to/rc/file.toml +``` +You can find an exemplary file in `rc/default_simulation.toml`. +The simulations will be saved as `.h5` files. \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 4f4c76d..4ef36a6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,7 @@ classifiers=[ dependencies = [ "astropy", "click", + "scikit-image", "h5py", "matplotlib", "numpy", diff --git a/radiosim/gauss.py b/radiosim/gauss.py index fd72f37..f3edb12 100644 --- a/radiosim/gauss.py +++ b/radiosim/gauss.py @@ -1,6 +1,8 @@ import numpy as np +from numpy.typing import ArrayLike from astropy.convolution import Gaussian2DKernel from astropy.nddata.utils import Cutout2D +from scipy.stats import norm, skewnorm def gauss(params: list, size: int): @@ -82,3 +84,60 @@ def rotgauss(x, y): return g return rotgauss(*np.indices((size, size))) + + +def skewed_gauss( + size: int, + x: float, + y: float, + amp: float, + width: float, + length: float, + a: float, +) -> ArrayLike: + """ + Generate a skewed 2d normal distribution. + + Parameters + ---------- + size: int + length of the square image + x: float + x coordinate of the center + y: float + y coordinate of the center + amp: float + maximal amplitude of the distribution + width: float + width of the distribution (perpendicular to the skewed function) + length: float + length of the distribution + a: float + skewness parameter + + Returns + ------- + skew_gauss: ArrayLike + two dimensional skewed gaussian distribution + """ + + def skewfunc(x: float, **args) -> ArrayLike: + func = skewnorm.pdf(x, **args) + func /= func.max() + return func + + vals = np.arange(size) + normal = norm.pdf(vals, loc=size / 2, scale=width).reshape(-1, 1) + normal /= normal.max() + + skew = skewfunc(vals, a=a, loc=size / 2, scale=length) + + skew_gauss = normal * skew + skew_gauss *= amp + skew_gauss[np.isclose(skew_gauss, 0)] = 0 + + skew_gauss = np.roll( + skew_gauss, (y - int(size / 2), x - int(size / 2)), axis=(0, 1) + ) + + return skew_gauss diff --git a/radiosim/mojave.py b/radiosim/mojave.py new file mode 100644 index 0000000..a9f43ab --- /dev/null +++ b/radiosim/mojave.py @@ -0,0 +1,601 @@ +import numpy as np +from numpy.typing import ArrayLike +from radiosim.gauss import skewed_gauss, twodgaussian +from radiosim.utils import _gen_date, _gen_vlba_obs_position +from scipy.stats import skewnorm, expon +from skimage.transform import swirl, rotate +from joblib import Parallel, delayed + + +def create_mojave(conf, rng): + """ + Create MOJAVE like sources. + + Parameters: + ----------- + conf : + loaded conf file + + Returns: + -------- + glx : ArrayLike + generated sources + """ + + threads = conf["threads"] + if threads == "none": + threads = 1 + elif not isinstance(threads, int): + raise ValueError("threads has to be int >0 or \"none\"") + elif threads <= 0: + raise ValueError("threads has to be int >0 or \"none\"") + + size = conf["img_size"] + + # calculate amount of each class to generate per bundle + bundle_size = conf["bundle_size"] + ratio = np.array(conf["class_ratio"]) + + # amount = np.array([n_comact, n_one_jet, n_two_jet]]) + amount = np.array((ratio / ratio.sum()) * bundle_size).astype(int) + while amount.sum() < bundle_size: + amount[amount.argmax()] += 1 + + jets = [] + compact = Parallel(n_jobs=threads)(delayed(gen_compact)(rng=child, size=size) for child in rng.spawn(amount[0])) + one_jet = Parallel(n_jobs=threads)(delayed(gen_one_jet)(rng=child, size=size) for child in rng.spawn(amount[1])) + two_jet = Parallel(n_jobs=threads)(delayed(gen_two_jet)(rng=child, size=size) for child in rng.spawn(amount[2])) + jets = np.array([*compact, *one_jet, *two_jet]) + + glx_class = [] + for glx_type, amt in zip([0, 1, 2], amount): + glx_class.extend([glx_type] * amt) + glx_class = np.array(glx_class) + + shuffler = np.arange(bundle_size) + rng.shuffle(shuffler) + + ra, dec = _gen_vlba_obs_position(rng=rng, size=bundle_size) + obs_dates = _gen_date(rng=rng, start_date="1995-01-01", size=bundle_size) + + data = [jets[shuffler], glx_class[shuffler], ra[shuffler], dec[shuffler], obs_dates] + data_name = ["galaxies", "galaxy_classes", "RA", "DEC", "date"] + + return data, data_name + +def gen_jet(size: int, amplitude: float, width: float, length: float, a: float) -> ArrayLike: + """ + Generate jet from skewed 2d normal distributiuon. + + Parameters + ---------- + size : int + length of the square image + x : float + x coordinate of the center + y : float + y coordinate of the center + amp : float + maximal amplitude of the distribution + width : float + width of the distribution (perpendicular to the skewed function) + length : float + length of the distribution + a : float + skewness parameter + + Returns + ------- + jet : ArrayLike + jet for source + """ + + jet = skewed_gauss(size=size, x=size/2, y=size/2, amp=amplitude, width=width, length=length, a=a) + + jet[np.isclose(jet, 0)] = 0 + + return jet + +def gen_shocks(jet: ArrayLike, shock_comp: int, length: float, width: float, sx: float, rng: "np.Generator") -> tuple[ArrayLike, tuple, int]: + """ + Generate equally spaced shocks in jet. + + Parameters + ---------- + jet : ArrayLike + generated jet without shock + shock_comp : int + max amount of shocks + length : float + length of jet + width : float + width of jet + sx : float + sx component of main gaussian + rng : np.Generator + numpy random generator + + Returns + ------- + jet : ArrayLike + jet with shock components + printout_information : tuple + information about each component + skipped : int + amount of skipped shock components + """ + ps_orientation = [] + ps_amp = [] + ps_x = [] + ps_y = [] + ps_sx = [] + ps_sy = [] + ps_rot = [] + ps_dist = [] + ps_a = [] + + mask = np.copy(jet) + mask[np.isclose(mask, 0)] = 0 + mask[mask > 0] = 1 + mask = mask.astype(bool) + + size = len(jet) + + s_length = 2 * length - sx + s_dist = s_length / shock_comp + s_dist0 = s_dist + skipped = 0 + for i in range(shock_comp): + # skip randomly + if rng.uniform(0, 1) > 2 / 3: + skipped += 1 + continue + + s_x = size / 2 + s_dist + s_y = size / 2 + s_amp = jet[int(s_y), int(s_x)] * rng.uniform(0.03, 0.6) + + s_sx = rng.uniform(*np.sort([2, width])) + s_sy = rng.uniform(*np.sort([s_sx, length / shock_comp])) + + s_rot = 0 + + s_a = rng.uniform(5, 9) + + jet += skewed_gauss(size, int(s_x), int(s_y), s_amp, s_sx, s_sy, s_a) + + ps_amp.append(s_amp) + ps_x.append(int(s_x)) + ps_y.append(int(s_y)) + ps_sx.append(s_sx) + ps_sy.append(s_sy) + ps_rot.append(s_rot) + ps_dist.append(s_dist) + ps_a.append(s_a) + ps_orientation.append("r") + + s_dist += s_dist0 + return ( + jet, + (ps_orientation, ps_amp, ps_x, ps_y, ps_sx, ps_sy, ps_rot, ps_dist, ps_a), + skipped, + ) + +def add_swirl(jet: ArrayLike, rng: "np.Generator", first_jet_params: tuple | None = None) -> tuple[ArrayLike, tuple]: + """ + Add swirl distortion to the jet. + + Parameters + ---------- + jet : ArrayLike + generated jet + rng : np.Generator + numpy random generator + first_jet_params : tuple | None, default: None + Use when generate two jet sources. + Tuple of parameters to generate similar swirl for second jet. + Contains the returned parameters from the first jet. + + Returns + ------- + swirled_jet : ArrayLike + swirl distorted input jet + parameters : tuple + parameters of the aplied swirl distortion + """ + size = len(jet) + if not first_jet_params: + strength = rng.normal(loc=0, scale=0.2) + radius = rng.uniform(50, 200) + center = np.array([size / 2, size / 2]) + else: + deviation = 1 + rng.normal(loc=0, scale=0.1, size=2) + strength = first_jet_params[0] * deviation[0] + radius = first_jet_params[1] * deviation[1] + center = first_jet_params[2] + + swirled_jet = swirl(jet, center=center, strength=strength, radius=radius) + + return swirled_jet, (strength, radius, center) + +def gen_two_jet(rng : int, size : int = 1024, printout : bool = False) -> ArrayLike: + """ + Generate a two jet galaxy. + + Parameters + ---------- + rng : np.Generator | int + numpy random generator or seed to generate random number generator + size : int, default: 1024 + size of the galaxy to be genrated + printout : bool, default: False + print out log information + + Returns + ------- + glx : ArrayLike + simulated two jet source + """ + if isinstance(rng, int): + rng = np.random.default_rng(seed=rng) + + amp_params = (1.715e-4, 1.985e-2) + width_params = (4.37, 22.71, 16.85) + length_params = (8.62, 69.60, 54.73) + left_deviation = 1 + rng.normal( + loc=0, scale=0.1, size=5 + ) # amp, length, width, strength, radius + + amp = expon.rvs(*amp_params, random_state=rng) + + r_width = skewnorm.rvs(*width_params, random_state=rng) / 10 # 10 + r_length = skewnorm.rvs(*length_params, random_state=rng) / 4 # 10 + r_jet_amp = amp * rng.power(0.7) + + l_width = r_width * left_deviation[2] + l_length = r_length * left_deviation[1] + l_jet_amp = r_jet_amp * left_deviation[0] + + dimensions = np.sort([r_width / 2, r_length / 5, l_width / 2, l_length / 5]) + + sx = rng.uniform(dimensions.min(), dimensions.max()) * 1.5 # too large + sy = rng.uniform(dimensions.min(), dimensions.max()) * 1.5 + while np.any( + [sx / sy < 1 / 2, sx / sy > 2] + ): # redraw if gaussian is too elliptical + sy = rng.uniform(dimensions.min(), dimensions.max()) * 1.5 + rot = rng.uniform(0, 2 * np.pi) + + a = 4 + + r_jet = gen_jet(size, r_jet_amp, r_width, r_length, a) + l_jet = gen_jet(size, l_jet_amp, l_width, l_length, a) + + shock_comp_r = int(np.round(rng.power(0.5), decimals=1) * 10) + shock_comp_l = int(np.round(rng.power(0.5), decimals=1) * 10) + + if shock_comp_r > 0: + r_jet, r_printout, r_skipped = gen_shocks(r_jet, shock_comp_r, r_length, r_width, sx, rng) + + if shock_comp_l > 0: + l_jet, l_printout, l_skipped = gen_shocks(l_jet, shock_comp_l, l_length, l_width, sx, rng) + + r_jet, r_swirl_comp = add_swirl(r_jet, rng=rng) + l_jet, l_swirl_comp = add_swirl(l_jet, rng=rng, first_jet_params=r_swirl_comp) + + glx = r_jet + np.flip(l_jet) + x = size / 2 + y = size / 2 + + gauss = twodgaussian([amp, size / 2, size / 2, sx, sy, rot], size) + glx += gauss + + angle = rng.uniform(0, 360) + + glx = rotate(glx, angle=angle) + + glx[np.isclose(glx, 0)] = 0 + + if printout: + print("left Jet:") + print( + f"width = {l_width:.3f}, length = {l_length:.3f}, \ + jet_amp = {l_jet_amp:.3f}, a = {a:.3f}" + ) + print("right Jet:") + print( + f"width = {r_width:.3f}, length = {r_length:.3f}, \ + jet_amp = {r_jet_amp:.3f}, a = {a:.3f}" + ) + print("Source:") + print( + f"amp = {amp:.3f}, x = {x:.0f}, y = {y:.0f}, sx = {sx:.3f}, \ + sy = {sy:.3f}, rot = {rot:.3f}" + ) + if shock_comp_r > 0: + print( + f"Right Shock, {shock_comp_r} Components, \ + {r_skipped} skipped:" + ) + for orient, amp, x, y, sx, sy, rot, dist, a in zip(*r_printout): + print( + f"amp = {amp:.6f}, x = {x:.0f}, y = {y:.0f}, \ + sx = {sx:.3f}, sy = {sy:.3f}, rot = {rot:.3f}, \ + dist = {dist:.3f}" + ) + if shock_comp_l > 0: + print( + f"Left Shock, {shock_comp_l} Components, \ + {l_skipped} skipped:" + ) + for orient, amp, x, y, sx, sy, rot, dist, a in zip(*l_printout): + print( + f"amp = {amp:.6f}, x = {x:.0f}, y = {y:.0f}, \ + sx = {sx:.3f}, sy = {sy:.3f}, rot = {rot:.3f}, \ + dist = {dist:.3f}" + ) + print("Swirl:") + print( + f"Right: strength = {r_swirl_comp[0]:.3f}, \ + radius = {r_swirl_comp[1]:.3f}, center = {r_swirl_comp[2]}" + ) + print( + f"Left: strength = {l_swirl_comp[0]:.3f}, \ + radius = {l_swirl_comp[1]:.3f}, center = {l_swirl_comp[2]}" + ) + + return glx + +def gen_one_jet(rng : int, size : int = 1024, printout : bool = False) -> ArrayLike: + """ + Generate a one jet galaxy. + + Parameters + ---------- + rng : np.Generator | int + numpy generator or seed to generate random number generator + size : int, default: 1024 + size of the galaxy to be genrated + printout : bool, default: False + print out debug information + + Returns + ------- + glx : ndarray + simulated one jet source + """ + if isinstance(rng, int): + rng = np.random.default_rng(seed=rng) + + amp_params = (1.715e-4, 1.985e-2) + width_params = (4.37, 22.71, 16.85) + length_params = (8.62, 69.60, 54.73) + + amp = expon.rvs(*amp_params, random_state=rng) + + width = skewnorm.rvs(*width_params, random_state=rng) / 12 # 10 + length = skewnorm.rvs(*length_params, random_state=rng) / 6.5 # 10 + jet_amp = amp * rng.power(0.7) + + dimensions = np.sort([width / 2, length / 5]) + + sx = rng.uniform(*dimensions) * 1.5 # too large + sy = rng.uniform(*dimensions) * 1.5 + while np.any( + [sx / sy < 1 / 2, sx / sy > 2] + ): # redraw if gaussian is too elliptical + sy = rng.uniform(dimensions.min(), dimensions.max()) * 1.5 + rot = rng.uniform(0, 2 * np.pi) + + a = 4 + + jet = gen_jet(size, jet_amp, width, length, a) + + shock_comp = int(np.round(rng.power(0.5), decimals=1) * 10) + + if shock_comp > 0: + jet, shock_printout, skipped = gen_shocks( + jet, shock_comp, length, width, sx, rng + ) + + jet, swirl_comp = add_swirl(jet, rng=rng) + + x = size / 2 + y = size / 2 + + gauss = twodgaussian([amp, size / 2, size / 2, sx, sy, rot], size) + glx = jet + gauss + + angle = rng.uniform(0, 360) + + glx = rotate(glx, angle=angle) + + glx[np.isclose(glx, 0)] = 0 + + if printout: + print("Skewed dist:") + print( + f"width = {width:.3f}, length = {length:.3f}, \ + jet_amp = {jet_amp:.3f}, a = {a:.3f}" + ) + print("Source:") + print( + f"amp = {amp:.3f}, x = {x:.0f}, y = {y:.0f}, sx = {sx:.3f}, \ + sy = {sy:.3f}, rot = {rot:.3f}" + ) + if shock_comp > 0: + print(f"{shock_comp} shock components, {skipped} skipped:") + for orient, amp, x, y, sx, sy, rot, dist, a in zip(shock_printout): + print( + f"amp = {amp:.6f}, x = {x:.0f}, y = {y:.0f}, \ + sx = {sx:.3f}, sy = {sy:.3f}, rot = {rot:.3f}, \ + dist = {dist:.3f}" + ) + print("Swirl:") + print( + f"strength = {swirl_comp[0]:.3f}, \ + radius = {swirl_comp[1]:.3f}, center = {swirl_comp[2]}" + ) + + return glx + +def gen_compact(rng : int, size : int = 1024, printout : bool = False) -> ArrayLike: + """Generate a compact galaxy. + + Parameters + ---------- + rng : np.Generator | int + numpy generator or seed to generate random number generator. + size : int, default: 1024 + size of the galaxy to be genrated + printout : bool, default: False + print out debug information. + + Returns + ------- + glx : ndarray + Simulated compact source. + """ + if isinstance(rng, int): + rng = np.random.default_rng(seed=rng) + + amp_params = (5.54e-05, 0.02) + width_params = (8.41, 3.04, 3.46) + length_params = (6.04, 3.02, 2.59) + + ncomps = rng.integers(1, 4) # draw amount of components + outflow = rng.uniform() > 0.5 # draw if outflow of the jet is generated + + amp = expon.rvs(*amp_params, random_state=rng) + + jet_amp = amp * rng.normal(2/3, 0.05) + + y = np.arange(size) + + width = skewnorm.rvs(*width_params, random_state=rng) * 0.8 + length = skewnorm.rvs(*length_params, random_state=rng) * 1 + + if outflow: + a = 4 + glx = gen_jet(size, jet_amp, width, length, a) + glx, swirl_comp = add_swirl(glx, rng=rng) + else: + glx = np.zeros((size, size)) + + dimensions = np.sort([width / 2, length / 4]) + dimensions[0] *= 1.2 + dimensions = np.sort(dimensions) + + i = 0 + comp_amp = [] + comp_x = [] + comp_y = [] + comp_sx = [] + comp_sy = [] + comp_rot = [] + vertical = True + while i < ncomps: + deviation = 1 + rng.normal(loc=0, scale=0.2, size=2) + if i != 0: + x = rng.uniform(-5, 5) + y = rng.uniform(-5, 5) + else: + x = 0 + y = 0 + if i == 0: + sx = rng.uniform(*dimensions) * 1.7 # too large + j = 0 + while sx > 15: + if j > 10: + sx /= 5 + break + sx = rng.uniform(*dimensions) * 1.7 + j += 1 + + sy = rng.uniform(*dimensions) * 1.7 + + j = 0 + while sy > 15: + if j > 10: + sy /= 5 + break + sy = rng.uniform(*dimensions) * 1.7 + j += 1 + + j = 0 + while np.any([sx/sy < 0.33, sx/sy>3]) or np.all([sx/sy > 0.75, sx/sy < 1.33]): + if j > 10: + sy = sx * rng.choice([rng.uniform(0.4, 0.8), rng.uniform(1.2, 2.5)]) + if sy > 15: + sy /= 5 + break + sy = rng.uniform(*dimensions) * 1.7 + k = 0 + while sy > 15: + if k > 10: + sy /= 5 + break + sy = rng.uniform(*dimensions) * 1.7 + k += 1 + j+=1 + else: + sx = comp_sx[-1] * deviation[0] + sy = comp_sy[-1] * deviation[1] + + if i == 0: + if sx > sy: + vertical = True + else: + vertical = False + + try: + rot = rng.uniform(comp_rot[-1] - 0.1, comp_rot[-1] + 0.1) + except IndexError: + if vertical: + rot = rng.uniform(-np.pi/4, np.pi/4) + else: + rot = rng.uniform(-np.pi/4+np.pi/2, np.pi/4+np.pi/2) + + c_amp = amp * rng.uniform(0.5, 1.2) + + gauss = twodgaussian([c_amp, size / 2 - x, size / 2 - y, sx, sy, rot], size) + glx += gauss + + comp_amp.append(c_amp) + comp_x.append(size / 2 - x) + comp_y.append(size / 2 - y) + comp_sx.append(sx) + comp_sy.append(sy) + comp_rot.append(rot) + + i += 1 + + angle = rng.uniform(0, 360) + glx = rotate(glx, angle=angle) + + shift = size // 2 - np.argwhere(glx == glx.max())[0] + + glx = np.roll(glx, shift=shift, axis=[0, 1]) + + glx[np.isclose(glx, 0)] = 0 + + glx /= glx.max() + glx *= amp + + if printout: + if outflow: + print("Skewed dist:") + print( + f"width = {width:.3f}, length = {length:.3f}, jet_amp = {jet_amp:.3f}, a = {a:.3f}" + ) + print("Swirl:") + print( + f"strength = {swirl_comp[0]:.3f}, radius = {swirl_comp[1]:.3f}, center = {swirl_comp[2]}" + ) + print("Sources:") + for cx, cy, csx, csy, cr in zip(comp_x, comp_y, comp_sx, comp_sy, comp_rot): + print( + f"amp = {amp:.3f}, x = {cx:.0f}, y = {cy:.0f}, sx = {csx:.3f}, sy = {csy:.3f}, rot = {cr:.3f}" + ) + + return glx diff --git a/radiosim/scripts/start_simulation.py b/radiosim/scripts/start_simulation.py index efa0158..0e56c5e 100644 --- a/radiosim/scripts/start_simulation.py +++ b/radiosim/scripts/start_simulation.py @@ -18,7 +18,7 @@ ), default="simulate", ) -def main(configuration_path, mode): +def main(configuration_path, mode) -> None: config = toml.load(configuration_path) conf = read_config(config) print(conf, "\n") diff --git a/radiosim/simulations.py b/radiosim/simulations.py index 6a662fb..8a401d6 100644 --- a/radiosim/simulations.py +++ b/radiosim/simulations.py @@ -5,9 +5,12 @@ add_noise, adjust_outpath, save_sky_distribution_bundle, + _save_mojave_bundle, ) from radiosim.jet import create_jet from radiosim.survey import create_survey +from radiosim.mojave import create_mojave +from numpy.random import default_rng def simulate_sky_distributions(conf): @@ -18,13 +21,27 @@ def simulate_sky_distributions(conf): ) -def create_sky_distribution(conf, opt: str): +def create_sky_distribution(conf, opt: str) -> None: + if conf["mode"] == "mojave": + seed = conf["seed"] + if seed == "none": + rng = default_rng() + elif isinstance(seed, int): + rng = default_rng(seed) + else: + raise TypeError("seed has to be int or \"none\"") + for _ in tqdm(range(conf["bundles_" + opt])): + path = adjust_outpath(conf["outpath"], "/samp_" + 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) + elif conf["mode"] == "mojave": + data, data_name = create_mojave(conf, rng) + _save_mojave_bundle(path, data=data, data_name=data_name) + continue else: click.echo("Given mode not found. Choose 'survey' or 'jet' in config file") @@ -35,5 +52,4 @@ 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) diff --git a/radiosim/utils.py b/radiosim/utils.py index 9280edb..888803c 100644 --- a/radiosim/utils.py +++ b/radiosim/utils.py @@ -6,6 +6,7 @@ import toml import click import numpy as np +from numpy.typing import ArrayLike from pathlib import Path from scipy import signal from astropy.convolution import Gaussian2DKernel @@ -206,8 +207,9 @@ def read_config(config): unpacked configurations """ sim_conf = {} - sim_conf["mode"] = config["general"]["mode"] sim_conf["seed"] = config["general"]["seed"] + sim_conf["mode"] = config["general"]["mode"] + sim_conf["threads"] = config["general"]["threads"] sim_conf["outpath"] = config["paths"]["outpath"] sim_conf["training_type"] = config["jet"]["training_type"] sim_conf["num_jet_components"] = config["jet"]["num_jet_components"] @@ -215,6 +217,7 @@ def read_config(config): sim_conf["num_sources"] = config["survey"]["num_sources"] sim_conf["class_distribution"] = config["survey"]["class_distribution"] sim_conf["scale_sources"] = config["survey"]["scale_sources"] + sim_conf["class_ratio"] = config["mojave"]["class_ratio"] sim_conf["bundles_train"] = config["image_options"]["bundles_train"] sim_conf["bundles_valid"] = config["image_options"]["bundles_valid"] sim_conf["bundles_test"] = config["image_options"]["bundles_test"] @@ -326,6 +329,41 @@ def save_sky_distribution_bundle(path, x, y, name_x="x", name_y="y"): hf.create_dataset(name_y, data=y) hf.close() +def _save_mojave_bundle(path : str, data : tuple, data_name : tuple) -> None: + """ + Write MOJAVE simulations created in analysis to h5 file. + + Parameters + ---------- + path: str + path to save file + data : tuple + data to store, e.g. galaxies, classes, coordinates, … + data_name : tuple + name of the data columns + """ + assert len(data) == len(data_name) + with h5py.File(path, "w") as hf: + for dat, name in zip(data, data_name): + hf.create_dataset(name, data=dat) + hf.close() + +def _gen_vlba_obs_position(rng, size : int) -> tuple[float, float]: + ra = rng.uniform(0, 360, size=size) + dec = rng.uniform(-30, 87, size=size) + return ra, dec + +def _gen_date(rng, start_date : str, size : int) -> ArrayLike: + from datetime import datetime + from h5py import string_dtype + # define fixed string length for h5py + utf8_type = string_dtype("utf-8", 10) + + start_date = np.datetime64(start_date) + now = np.datetime64(datetime.now().date()) + delta = now - start_date + dates = start_date + rng.integers(delta.astype(int), size=size) + return dates.astype(utf8_type) def cart2pol(x: float, y: float): """ diff --git a/rc/default_simulation.toml b/rc/default_simulation.toml index 691343b..983ccac 100644 --- a/rc/default_simulation.toml +++ b/rc/default_simulation.toml @@ -5,7 +5,8 @@ title = "Simulation configuration" [general] quiet = true seed = "none" # "none" or int -mode = "jet" # jet or survey +mode = "jet" # jet or survey or mojave +threads = "none" # "none" or int > 0 [paths] outpath = "./build/example_data/" @@ -21,6 +22,9 @@ num_sources = 20 class_distribution = [2, 1, 2] # average share of each source in the images scale_sources = true +[mojave] +class_ratio = [5, 5, 1] # rel. amount of each class [compact, one_jet, two_jets] + [image_options] bundles_train = 1 bundles_valid = 1 diff --git a/setup.py b/setup.py deleted file mode 100644 index d5e43c3..0000000 --- a/setup.py +++ /dev/null @@ -1,45 +0,0 @@ -from setuptools import setup, find_packages - -setup( - name="radiosim", - 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", - author_email="kevin3.schmidt@tu-dortmund.de", - license="MIT", - packages=find_packages(), - install_requires=[ - "numpy", - "matplotlib", - "h5py", - "toml", - "click", - "tqdm", - "scipy", - "astropy", - "pathlib", - "pytest", - "opencv-python", - ], - setup_requires=["pytest-runner"], - tests_require=["pytest"], - zip_safe=False, - entry_points={ - "console_scripts": [ - "radiosim = radiosim.scripts.start_simulation:main", - ], - }, - classifiers=[ - "Development Status :: 3 - Alpha", - "Intended Audience :: Science/Research", - "License :: OSI Approved :: MIT License", - "Natural Language :: English", - "Operating System :: OS Independent", - "Programming Language :: Python", - "Programming Language :: Python :: 3.9", - "Programming Language :: Python :: 3 :: Only", - "Topic :: Scientific/Engineering :: Astronomy", - "Topic :: Scientific/Engineering :: Physics", - ], -)