From 3d0535ce04f352641a5acf1647440f61ded4367d Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Wed, 26 Jul 2023 11:15:21 +0200 Subject: [PATCH 01/63] Use nearest-neighbour interpolation in regions where extrapolation is required. (#285) * Fix version drift. * Bump to 0.2.0 * Use nearest-neighbour interpolation for points requiring extrapolation. --- pyproject.toml | 2 +- quartical/interpolation/interpolants.py | 23 ++++++++++++++++++++++- tbump.toml | 2 +- 3 files changed, 24 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index a2196b67..68c1f4ee 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "quartical" -version = "0.1.10" +version = "0.2.0" description = "Fast and flexible calibration suite for radio interferometer data." repository = "https://github.com/ratt-ru/QuartiCal" documentation = "https://quartical.readthedocs.io" diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 0c384645..07771157 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -35,11 +35,32 @@ def linear2d_interpolate_gains(source_xds, target_xds): if i_f_dim > 1: interp_axes[i_f_axis] = target_xds[t_f_axis].data - target_xda = source_xds.params.interp( + # NOTE: The below is the path of least resistance but may not be the most + # efficient method for this mixed-mode interpolation - it may be possible + # to do better using multiple RegularGridInterpoator objects. + + # Interpolate using linear interpolation, filling points outside the + # domain with NaNs. + in_domain_xda = source_xds.params.interp( + interp_axes, + kwargs={"fill_value": np.nan} + ) + + # Interpolate using nearest-neighbour interpolation, extrapolating points + # outside the domain. + out_domain_xda = source_xds.params.interp( interp_axes, + method="nearest", kwargs={"fill_value": "extrapolate"} ) + # Combine the linear and nearest neighbour interpolation done above i.e. + # use use linear interpolation inside the domain and nearest-neighbour + # interpolation anywhere extrapolation was required. + target_xda = in_domain_xda.where( + da.isfinite(in_domain_xda), out_domain_xda + ) + if i_t_dim == 1: target_xda = target_xda.reindex( {i_t_axis: target_xds[t_t_axis].data}, diff --git a/tbump.toml b/tbump.toml index e481d3d3..78ae3488 100644 --- a/tbump.toml +++ b/tbump.toml @@ -2,7 +2,7 @@ github_url = "https://github.com/ratt-ru/QuartiCal/" [version] -current = "0.1.11" +current = "0.2.0" # Example of a semver regexp. # Make sure this matches current_version before From 6b22750e309dc3d715c2512f3afc38218809dc66 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 28 Jul 2023 12:23:54 +0200 Subject: [PATCH 02/63] Add interface for degrid models. --- quartical/config/preprocess.py | 49 +++++++++++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/quartical/config/preprocess.py b/quartical/config/preprocess.py index d7fa2081..ce75f009 100644 --- a/quartical/config/preprocess.py +++ b/quartical/config/preprocess.py @@ -6,15 +6,31 @@ import os.path from dataclasses import dataclass from typing import List, Dict, Set, Any +from ast import literal_eval sky_model_nt = namedtuple("sky_model_nt", ("name", "tags")) +degrid_model_nt = namedtuple( + "degrid_model_nt", + ( + "name", + "nxo", + "nyo", + "cellxo", + "cellyo", + "x0o", + "y0o", + "ipi", + "cpi" + ) +) @dataclass class Ingredients: model_columns: Set[Any] sky_models: Set[sky_model_nt] + degrid_models: Set[degrid_model_nt] @dataclass @@ -49,6 +65,8 @@ def transcribe_recipe(user_recipe): model_columns = set() sky_models = set() + degrid_models = set() + instructions = {} # Strip accidental whitepsace from input recipe and splits on ":". @@ -92,6 +110,18 @@ def transcribe_recipe(user_recipe): sky_models.add(sky_model) instructions[recipe_index].append(sky_model) + elif ".mds" in ingredient: + + filename, _, options = ingredient.partition("@") + options = literal_eval(options) # Add fail on missing option. + + if not os.path.exists(filename): + raise FileNotFoundError("{} not found.".format(filename)) + + degrid_model = degrid_model_nt(filename, *options) + degrid_models.add(degrid_model) + instructions[recipe_index].append(degrid_model) + elif ingredient != "": model_columns.add(ingredient) instructions[recipe_index].append(ingredient) @@ -99,19 +129,32 @@ def transcribe_recipe(user_recipe): else: instructions[recipe_index].append(ingredient) + # TODO: Add message to log. logger.info("The following model sources were obtained from " "--input-model-recipe: \n" " Columns: {} \n" - " Sky Models: {}", + " Sky Models: {} \n" + " Degrid Models: {}", model_columns or 'None', - {sm.name for sm in sky_models} or 'None') + {sm.name for sm in sky_models} or 'None', + {dm.name for dm in degrid_models} or 'None') # Generate a named tuple containing all the information required to # build the model visibilities. - model_recipe = Recipe(Ingredients(model_columns, sky_models), instructions) + model_recipe = Recipe( + Ingredients( + model_columns, + sky_models, + degrid_models + ), + instructions + ) if model_recipe.ingredients.sky_models: logger.info("Recipe contains sky models - enabling prediction step.") + if model_recipe.ingredients.degrid_models: + logger.info("Recipe contains degrid models - enabling degridding.") + return model_recipe From a6b6accd934bbe1c53044f30444ff4cbf660c65b Mon Sep 17 00:00:00 2001 From: landmanbester Date: Fri, 28 Jul 2023 16:16:47 +0200 Subject: [PATCH 03/63] Add optional label and single field selection to backup app --- quartical/apps/backup.py | 29 ++++++++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/quartical/apps/backup.py b/quartical/apps/backup.py index fcb44a69..58e65fce 100644 --- a/quartical/apps/backup.py +++ b/quartical/apps/backup.py @@ -10,7 +10,7 @@ def backup(): parser = argparse.ArgumentParser( description='Backup any Measurement Set column to zarr. Backups will ' - 'be labelled automatically using the current datetime, ' + 'be labelled the passed in label (defaults to datetime), ' 'the Measurement Set name and the column name.' ) @@ -33,19 +33,32 @@ def backup(): type=str, help='Name of column to be backed up.' ) + parser.add_argument( + 'label', + type=str, + help='A label for the backup. Default to datetime.' + ) parser.add_argument( '--nthread', type=int, default=1, help='Number of threads to use.' ) + parser.add_argument( + '--fieldid', + type=int, + help='Field ID to back up.' + ) args = parser.parse_args() ms_name = args.ms_path.full_path.rsplit("/", 1)[1] column_name = args.column_name - timestamp = time.strftime("%Y%m%d-%H%M%S") + if args.label is not None: + label = args.label + else: + label = time.strftime("%Y%m%d-%H%M%S") # This call exists purely to get the relevant shape and dtype info. data_xds_list = xds_from_storage_ms( @@ -55,6 +68,11 @@ def backup(): group_cols=("FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER"), ) + if args.fieldid is not None: + for ds in data_xds_list: + if ds.FIELD_ID != args.fieldid: + data_xds_list.pop(ds) + # Compute appropriate chunks (256MB by default) to keep zarr happy. chunks = [chunk_by_size(xds[column_name]) for xds in data_xds_list] @@ -67,9 +85,14 @@ def backup(): chunks=chunks ) + if args.fieldid is not None: + for ds in data_xds_list: + if ds.FIELD_ID != args.fieldid: + data_xds_list.pop(ds) + bkp_xds_list = xds_to_zarr( data_xds_list, - f"{args.zarr_dir.url}::{timestamp}-{ms_name}-{column_name}.bkp.qc", + f"{args.zarr_dir.url}::{label}-{ms_name}-{column_name}.bkp.qc", ) dask.compute(bkp_xds_list, num_workers=args.nthread) From e1e3dcb48ccbb8933beee66140a6bd9c75d4684d Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 28 Jul 2023 17:40:33 +0200 Subject: [PATCH 04/63] Commit initial, semi-working code for generating model visibilities from pfb-style models. --- quartical/config/preprocess.py | 2 +- quartical/data_handling/degridder.py | 196 +++++++++++++++++++++++ quartical/data_handling/model_handler.py | 40 ++++- 3 files changed, 230 insertions(+), 8 deletions(-) create mode 100644 quartical/data_handling/degridder.py diff --git a/quartical/config/preprocess.py b/quartical/config/preprocess.py index ce75f009..ea11d82b 100644 --- a/quartical/config/preprocess.py +++ b/quartical/config/preprocess.py @@ -155,6 +155,6 @@ def transcribe_recipe(user_recipe): logger.info("Recipe contains sky models - enabling prediction step.") if model_recipe.ingredients.degrid_models: - logger.info("Recipe contains degrid models - enabling degridding.") + logger.info("Recipe contains degrid models - enabling degridding.") return model_recipe diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py new file mode 100644 index 00000000..24711897 --- /dev/null +++ b/quartical/data_handling/degridder.py @@ -0,0 +1,196 @@ +from collections import defaultdict +import numpy as np +import dask.array as da +from daskms.experimental.zarr import xds_from_zarr +from scipy.interpolate import RegularGridInterpolator +import sympy as sm +from sympy.utilities.lambdify import lambdify +from sympy.parsing.sympy_parser import parse_expr +from ducc0.wgridder.experimental import dirty2vis +from quartical.utils.collections import freeze_default_dict + + +def _degrid(time, freq, uvw, model): + + name, nxo, nyo, cellxo, cellyo, x0o, y0o, ipi, cpi = model + + # TODO: Dodgy pattern, as this will produce dask arrays which will + # need to be evaluated i.e. task-in-task which is a bad idea. OK for + # the purposes of a prototype. + model_xds = xds_from_zarr(model.name)[0] + + # TODO: We want to go from the model to an image cube appropriate for this + # chunks of data. The dimensions of the image cube will be determined + # by the time and freqency dimensions of the chunk in conjunction with + # the integrations per chunk and channels per chunk arguments. + + utime, utime_inv = np.unique(time, return_inverse=True) + n_utime = utime.size + ipi = ipi or n_utime # Catch zero case. NOTE: -1? + n_mean_times = int(np.ceil(n_utime / ipi)) + + n_freq = freq.size + n_mean_freqs = int(np.ceil(n_freq / cpi)) + + # Let's start off by simply reconstructing the model image as generated + # by pfb-clean. + + nxi, nyi = model_xds.npix_x, model_xds.npix_y + image_in = np.zeros((nxi, nyi), dtype=float) + params = sm.symbols(('t', 'f')) + params += sm.symbols(tuple(model_xds.params.values)) + symexpr = parse_expr(model_xds.parametrisation) + modelf = lambdify(params, symexpr) + texpr = parse_expr(model_xds.texpr) + tfunc = lambdify(params[0], texpr) + fexpr = parse_expr(model_xds.fexpr) + ffunc = lambdify(params[1], fexpr) + Ix = model_xds.location_x.values + Iy = model_xds.location_y.values + coeffs = model_xds.coefficients.values + + vis = np.empty((time.size, freq.size, 4), dtype=np.complex128) + + for ti in range(n_mean_times): + for fi in range(n_mean_freqs): + + freq_sel = slice(fi * cpi, (fi + 1) * cpi) + degrid_freq = freq[freq_sel] + mean_freq = degrid_freq.mean() + + time_sel = slice(ti * ipi, (ti + 1) * ipi) + degrid_time = utime[time_sel] + mean_time = degrid_time.mean() + + image_in[Ix, Iy] = modelf( + tfunc(mean_time), ffunc(mean_freq), *coeffs + ) + + # NOTE: Select out appropriate rows for ipi and make selection + # consistent. + row_sel = slice(None) + + # Degrid + dirty2vis( + vis=vis[row_sel, freq_sel, 0], + uvw=uvw, + freq=degrid_freq, + dirty=image_in, + pixsize_x=model_xds.cell_rad_x, # Should be output value. + pixsize_y=model_xds.cell_rad_y, # Should be output value. + center_x=model_xds.center_x, # Should be output value. + center_y=model_xds.center_y, # Should be output value. + epsilon=1e-7, # TODO: Is this too high? + do_wgridding=True, # Should be ok to leave True. + divide_by_n=False, # Until otherwise informed. + nthreads=6 # Should be equivalent to solver threads. + ) + + # Zero the image array between image slices. + image_in[:, :] = 0 + + # Degridder only produces I - will need to be more sophisticated. + vis[..., -1] = vis[..., 0] + + return vis + + # TODO: This was omitted for the sake of simplicity but we ultimately will + # want this functionality. Need to be very cautious with regard to which + # parameters get used. + + cellxi, cellyi = model_xds.cell_rad_x, model_xds.cell_rad_y + x0i, y0i = model_xds.center_x, model_xds.center_y + + xin = (-(nxi//2) + np.arange(nxi))*cellxi + x0i + yin = (-(nyi//2) + np.arange(nyi))*cellyi + y0i + xo = (-(nxo//2) + np.arange(nxo))*cellxo + x0o + yo = (-(nyo//2) + np.arange(nyo))*cellyo + y0o + + # how many pixels to pad by to extrapolate with zeros + xldiff = xin.min() - xo.min() + if xldiff > 0.0: + npadxl = int(np.ceil(xldiff/cellxi)) + else: + npadxl = 0 + yldiff = yin.min() - yo.min() + if yldiff > 0.0: + npadyl = int(np.ceil(yldiff/cellyi)) + else: + npadyl = 0 + + xudiff = xo.max() - xin.max() + if xudiff > 0.0: + npadxu = int(np.ceil(xudiff/cellxi)) + else: + npadxu = 0 + yudiff = yo.max() - yin.max() + if yudiff > 0.0: + npadyu = int(np.ceil(yudiff/cellyi)) + else: + npadyu = 0 + + do_pad = npadxl > 0 + do_pad |= npadxu > 0 + do_pad |= npadyl > 0 + do_pad |= npadyu > 0 + if do_pad: + image_in = np.pad( + image_in, + ((npadxl, npadxu), (npadyl, npadyu)), + mode='constant' + ) + + xin = (-(nxi//2+npadxl) + np.arange(nxi + npadxl + npadxu))*cellxi + x0i + nxi = nxi + npadxl + npadxu + yin = (-(nyi//2+npadyl) + np.arange(nyi + npadyl + npadyu))*cellyi + y0i + nyi = nyi + npadyl + npadyu + + do_interp = cellxi != cellxo + do_interp |= cellyi != cellyo + do_interp |= x0i != x0o + do_interp |= y0i != y0o + do_interp |= nxi != nxo + do_interp |= nyi != nyo + if do_interp: + interpo = RegularGridInterpolator((xin, yin), image_in, + bounds_error=True, method='linear') + xx, yy = np.meshgrid(xo, yo, indexing='ij') + return interpo((xx, yy)) + # elif (nxi != nxo) or (nyi != nyo): + # # only need the overlap in this case + # _, idx0, idx1 = np.intersect1d(xin, xo, assume_unique=True, return_indices=True) + # _, idy0, idy1 = np.intersect1d(yin, yo, assume_unique=True, return_indices=True) + # return image[idx0, idy0] + else: + return image_in + + +def degrid(data_xds_list, model_vis_recipe, ms_path, model_opts): + + degrid_models = model_vis_recipe.ingredients.degrid_models + + degrid_list = [] + + for data_xds in data_xds_list: + + model_vis = defaultdict(list) + + for degrid_model in degrid_models: + + degrid_vis = da.blockwise( + _degrid, ("rowlike", "chan", "corr"), + data_xds.TIME.data, ("rowlike",), + data_xds.CHAN_FREQ.data, ("chan",), + data_xds.UVW.data, ("rowlike", "uvw"), + degrid_model, None, + concatenate=True, + align_arrays=False, + meta=np.empty([0, 0, 0], dtype=np.complex128), + new_axes={"corr": 4} # TODO: Shouldn't be hardcoded. + ) + + model_vis[degrid_model].append(degrid_vis) + + degrid_list.append(freeze_default_dict(model_vis)) + + return degrid_list diff --git a/quartical/data_handling/model_handler.py b/quartical/data_handling/model_handler.py index 9871c5ad..e5b6943e 100644 --- a/quartical/data_handling/model_handler.py +++ b/quartical/data_handling/model_handler.py @@ -2,6 +2,7 @@ import dask.array as da import numpy as np from quartical.data_handling.predict import predict +from quartical.data_handling.degridder import degrid from quartical.data_handling.angles import apply_parangles from quartical.config.preprocess import IdentityRecipe, Ingredients from quartical.utils.array import flat_ident_like @@ -36,13 +37,38 @@ def add_model_graph( # Generates a predicition scheme (graph) per-xds. If no predict is # required, it is a list of empty dictionaries. - if model_vis_recipe.ingredients.sky_models: - predict_schemes = predict(data_xds_list, - model_vis_recipe, - ms_path, - model_opts) - else: - predict_schemes = [{}]*len(data_xds_list) + # TODO: Add handling for mds inputs. This will need to read the model + # and figure out the relevant steps to take. + + predict_required = bool(model_vis_recipe.ingredients.sky_models) + degrid_required = bool(model_vis_recipe.ingredients.degrid_models) + + # TODO: Ensure that things work correctly when we have a mixture of the + # below. + + predict_schemes = [{}]*len(data_xds_list) + + if predict_required: + rime_schemes = predict( + data_xds_list, + model_vis_recipe, + ms_path, + model_opts + ) + predict_schemes = [ + {**ps, **rs} for ps, rs in zip(predict_schemes, rime_schemes) + ] + + if degrid_required: + degrid_schemes = degrid( + data_xds_list, + model_vis_recipe, + ms_path, + model_opts + ) + predict_schemes = [ + {**ps, **ds} for ps, ds in zip(predict_schemes, degrid_schemes) + ] # Special case: in the event that we have an IdentityRecipe, modify the # datasets and model appropriately. From 40f5029e8339fb833ca07883065b44926560c5ad Mon Sep 17 00:00:00 2001 From: landmanbester Date: Sun, 30 Jul 2023 15:07:16 +0200 Subject: [PATCH 05/63] remove item instead of pop@index --- quartical/apps/backup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/quartical/apps/backup.py b/quartical/apps/backup.py index 58e65fce..ad3696e2 100644 --- a/quartical/apps/backup.py +++ b/quartical/apps/backup.py @@ -71,7 +71,7 @@ def backup(): if args.fieldid is not None: for ds in data_xds_list: if ds.FIELD_ID != args.fieldid: - data_xds_list.pop(ds) + data_xds_list.remove(ds) # Compute appropriate chunks (256MB by default) to keep zarr happy. chunks = [chunk_by_size(xds[column_name]) for xds in data_xds_list] @@ -88,7 +88,7 @@ def backup(): if args.fieldid is not None: for ds in data_xds_list: if ds.FIELD_ID != args.fieldid: - data_xds_list.pop(ds) + data_xds_list.remove(ds) bkp_xds_list = xds_to_zarr( data_xds_list, From d7e69e9821675dc056fbc8bcf69e227691baecaa Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 09:45:05 +0200 Subject: [PATCH 06/63] Add ducc0 and sympy as extras under degrid banner in pyproject.toml. --- pyproject.toml | 5 +++++ quartical/data_handling/model_handler.py | 12 +++++++++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 68c1f4ee..ceee48d2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,6 +40,11 @@ pytest = "^7.3.1" omegaconf = "^2.3.0" colorama = "^0.4.6" stimela = "2.0rc4" +ducc0 = "^0.31.0" +sympy = "^1.12" + +[tool.poetry.extras] +degrid = ["ducc0", "sympy"] [tool.poetry.scripts] goquartical = 'quartical.executor:execute' diff --git a/quartical/data_handling/model_handler.py b/quartical/data_handling/model_handler.py index e5b6943e..9c237809 100644 --- a/quartical/data_handling/model_handler.py +++ b/quartical/data_handling/model_handler.py @@ -2,7 +2,6 @@ import dask.array as da import numpy as np from quartical.data_handling.predict import predict -from quartical.data_handling.degridder import degrid from quartical.data_handling.angles import apply_parangles from quartical.config.preprocess import IdentityRecipe, Ingredients from quartical.utils.array import flat_ident_like @@ -60,6 +59,17 @@ def add_model_graph( ] if degrid_required: + + try: + from quartical.data_handling.degridder import degrid + except ImportError: + raise ImportError( + "QuartiCal was unable to import the degrid module. This may " + "indicate that QuartiCal was installed without the necessary " + "extras. Please try 'pip install quartical[degrid]'. If the " + "error persists, please raise an issue." + ) + degrid_schemes = degrid( data_xds_list, model_vis_recipe, From d8b6f3afbc5bf1c4ab99bf0a1b49769dfed4154c Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 31 Jul 2023 12:15:41 +0200 Subject: [PATCH 07/63] do not .remove() from xds_list --- quartical/apps/backup.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/quartical/apps/backup.py b/quartical/apps/backup.py index ad3696e2..3f36024f 100644 --- a/quartical/apps/backup.py +++ b/quartical/apps/backup.py @@ -68,13 +68,17 @@ def backup(): group_cols=("FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER"), ) + xdso = [] if args.fieldid is not None: for ds in data_xds_list: - if ds.FIELD_ID != args.fieldid: - data_xds_list.remove(ds) + if ds.FIELD_ID == args.fieldid: + # data_xds_list.remove(ds) # this fails with reindexing error + xdso.append(ds) + else: + xdso = data_xds_list # Compute appropriate chunks (256MB by default) to keep zarr happy. - chunks = [chunk_by_size(xds[column_name]) for xds in data_xds_list] + chunks = [chunk_by_size(xds[column_name]) for xds in xdso] # Repeat of above call but now with correct chunking information. data_xds_list = xds_from_storage_ms( @@ -85,13 +89,17 @@ def backup(): chunks=chunks ) + xdso = [] if args.fieldid is not None: for ds in data_xds_list: - if ds.FIELD_ID != args.fieldid: - data_xds_list.remove(ds) + if ds.FIELD_ID == args.fieldid: + # data_xds_list.remove(ds) + xdso.append(ds) + else: + xdso = data_xds_list bkp_xds_list = xds_to_zarr( - data_xds_list, + xdso, f"{args.zarr_dir.url}::{label}-{ms_name}-{column_name}.bkp.qc", ) From 0b26002afa922af766ed0d5d71a438a8798e3d9b Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 14:59:56 +0200 Subject: [PATCH 08/63] Add/change functionlity to enable new-style models. --- quartical/config/__init__.py | 22 ++- quartical/config/argument_schema.yaml | 6 + quartical/config/config_classes.py | 9 ++ quartical/config/external.py | 22 ++- quartical/config/helper.py | 11 +- quartical/config/internal.py | 10 ++ quartical/config/model_component_schema.yaml | 61 +++++++++ quartical/config/preprocess.py | 134 ++++++++++++++++++- quartical/executor.py | 11 +- 9 files changed, 274 insertions(+), 12 deletions(-) create mode 100644 quartical/config/model_component_schema.yaml diff --git a/quartical/config/__init__.py b/quartical/config/__init__.py index b4167112..b66a5315 100644 --- a/quartical/config/__init__.py +++ b/quartical/config/__init__.py @@ -27,8 +27,8 @@ class _GainSchema(object): gain: Dict[str, Parameter] - # The gain section is loaded explicitly, since we need to form up multiple - # instances. + # The gain and model sections are loaded explicitly, since we need to form + # up multiple instances. gain_schema = oc.merge( oc.structured(_GainSchema), oc.load(f"{dirname}/gain_schema.yaml") @@ -42,3 +42,21 @@ class _GainSchema(object): bases=(BaseConfigSection,), post_init=POST_INIT_MAP['gain'] ) + + @dataclass + class _ModelComponentSchema(object): + model_component: Dict[str, Parameter] + + model_component_schema = oc.merge( + oc.structured(_ModelComponentSchema), + oc.load(f"{dirname}/model_component_schema.yaml") + ) + model_component_schema = model_component_schema.model_component + + # Create model dataclass. + ModelComponent = schema_utils.schema_to_dataclass( + model_component_schema, + "ModelComponent", + bases=(BaseConfigSection,), + post_init=POST_INIT_MAP['model_component'] + ) diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index 5e5ad044..28644749 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -119,6 +119,12 @@ input_model: given by the dE tagged clusters in LSM1. Leaving this value unset (the default) will use an identity model. + advanced_recipe: + dtype: Optional[str] + info: + Optional advanced recipe specification for use with degridder/more + complicated modes of model construction. + beam: dtype: Optional[str] info: diff --git a/quartical/config/config_classes.py b/quartical/config/config_classes.py index 5f67e0f7..9baffab9 100644 --- a/quartical/config/config_classes.py +++ b/quartical/config/config_classes.py @@ -111,6 +111,14 @@ def __input_model_post_init__(self): self.__validate_choices__() self.__validate_element_choices__() + assert not (self.recipe and self.advanced_recipe), \ + "recipe and advanced_recipe are mutually exclusive." + + +def __model_component_post_init__(self): + self.__validate_choices__() + self.__validate_element_choices__() + def __output_post_init__(self): @@ -182,6 +190,7 @@ def __gain_post_init__(self): POST_INIT_MAP = { "input_ms": __input_ms_post_init__, "input_model": __input_model_post_init__, + "model_component": __model_component_post_init__, "output": __output_post_init__, "mad_flags": __mad_flags_post_init__, "solver": __solver_post_init__, diff --git a/quartical/config/external.py b/quartical/config/external.py index ce426e05..5c3de452 100644 --- a/quartical/config/external.py +++ b/quartical/config/external.py @@ -1,8 +1,9 @@ +import re from dataclasses import make_dataclass from omegaconf import OmegaConf as oc from typing import Dict, Any from scabha.cargo import Parameter -from quartical.config import Gain, BaseConfig, gain_schema +from quartical.config import Gain, ModelComponent, BaseConfig, gain_schema def finalize_structure(additional_config): @@ -18,9 +19,26 @@ def finalize_structure(additional_config): # Use the default terms if no alternative is specified. terms = terms or BaseConfig.solver.terms + recipe = None + models = [] # No components by default. + + # Get last specified version of input_model.recipe. + for cfg in additional_config[::-1]: + recipe = oc.select(cfg, "input_model.advanced_recipe") + if recipe is not None: + ingredients = re.split(r'([\+~:])', recipe) + ingredients = [ + i for i in ingredients if not bool(re.search(r'([\+~:])', i)) + ] + models = set(i.split("@")[0] for i in ingredients) + break + FinalConfig = make_dataclass( "FinalConfig", - [(t, Gain, Gain()) for t in terms], + [ + *[(m, ModelComponent, ModelComponent()) for m in models], + *[(t, Gain, Gain()) for t in terms] + ], bases=(BaseConfig,) ) diff --git a/quartical/config/helper.py b/quartical/config/helper.py index f4885da6..00ebecdb 100644 --- a/quartical/config/helper.py +++ b/quartical/config/helper.py @@ -79,8 +79,15 @@ def help(): if len(sys.argv) != 1 and not help_arg: return - # Include a generic gain term in the help config. - additional_config = [oc.from_dotlist(["solver.terms=['gain']"])] + # Include a generic gain term and model component in the help config. + additional_config = [ + oc.from_dotlist( + [ + "input_model.advanced_recipe=model_component", + "solver.terms=['gain']" + ] + ) + ] HelpConfig = finalize_structure(additional_config) if len(sys.argv) == 1 or help_arg == "help": diff --git a/quartical/config/internal.py b/quartical/config/internal.py index e86e6a9f..a8b06231 100644 --- a/quartical/config/internal.py +++ b/quartical/config/internal.py @@ -1,5 +1,15 @@ from daskms.fsspec_store import DaskMSStore from quartical.gains import TERM_TYPES +from quartical.config import ModelComponent + + +def get_component_dict(opts): + + return { + k: getattr(opts, k) + for k in opts.__dataclass_fields__.keys() + if isinstance(getattr(opts, k), ModelComponent) + } def gains_to_chain(opts): diff --git a/quartical/config/model_component_schema.yaml b/quartical/config/model_component_schema.yaml new file mode 100644 index 00000000..b1629c23 --- /dev/null +++ b/quartical/config/model_component_schema.yaml @@ -0,0 +1,61 @@ +model_component: + path_or_name: + dtype: str + required: true + info: + Path to model/name of column. + type: + dtype: str + choices: + - mds # Replace with awesome acronym eventually. + - tigger-lsm + - column + default: column + info: + Type of model component + tags: + dtype: Optional[List[str]] + info: + Tag for use in the predict. + region: + dtype: Optional[str] + info: + Name of region file to use in degridder. + npix_x: + dtype: Optional[int] + info: + Image x size in pixels for use in degridding. + npix_y: + dtype: Optional[int] + info: + Image y size in pixels for use in degridding. + cellsize_x: + dtype: Optional[float] + info: + Pixel x cellsize in radians. + cellsize_y: + dtype: Optional[float] + info: + Pixel y cellsize in radians. + centre_x: + dtype: int + default: 0 + info: + x coordinate of central pixel. + centre_y: + dtype: int + default: 0 + info: + y coordinate of central pixel. + integrations_per_image: + dtype: int + default: 0 + info: + Number of integrations per image to use in degridding. The default + (zero) is all times in a chunk. + channels_per_image: + dtype: int + default: 0 + info: + Number of channels per image to use in degridding. The default + (zero) is all channels in a chunk. diff --git a/quartical/config/preprocess.py b/quartical/config/preprocess.py index ea11d82b..c016067d 100644 --- a/quartical/config/preprocess.py +++ b/quartical/config/preprocess.py @@ -44,7 +44,7 @@ class IdentityRecipe(Recipe): pass -def transcribe_recipe(user_recipe): +def transcribe_legacy_recipe(user_recipe): """Interpret the model recipe string. Given the config object, create an internal recipe implementing the user @@ -72,10 +72,6 @@ def transcribe_recipe(user_recipe): # Strip accidental whitepsace from input recipe and splits on ":". input_recipes = user_recipe.replace(" ", "").split(":") - if input_recipes == ['']: - raise ValueError("No model recipe was specified. Please set/check " - "--input-model-recipe.") - for recipe_index, recipe in enumerate(input_recipes): instructions[recipe_index] = [] @@ -158,3 +154,131 @@ def transcribe_recipe(user_recipe): logger.info("Recipe contains degrid models - enabling degridding.") return model_recipe + + +def transcribe_recipe(user_recipe, model_components): + """Interpret the model recipe string. + + Given the config object, create an internal recipe implementing the user + specified recipe. + + Args: + model_opts: An ModelInputs configuration object. + + Returns: + model_Recipe: A Recipe object. + """ + + if user_recipe is None: + logger.warning( + "input_model.recipe was not supplied. Assuming identity model." + ) + return IdentityRecipe(Ingredients(set(), set()), dict()) + + model_columns = set() + sky_models = set() + degrid_models = set() + + instructions = {} + + # Strip accidental whitepspace from input recipe and splits on ":". + input_recipes = user_recipe.replace(" ", "").split(":") + + for recipe_index, recipe in enumerate(input_recipes): + + instructions[recipe_index] = [] + + # A raw string is required to avoid insane escape characters. Splits + # on understood operators, ~ for subtract, + for add. + + ingredients = re.split(r'([\+~])', recipe) + + # Behaviour of re.split guarantees every second term is either a column + # or .lsm. This may lead to the first element being an empty string. + + # Split the ingredients into operations and model sources. We preserve + # empty strings in the recipe to avoid more complicated code elsewhere. + + for ingredient in ingredients: + + if ingredient in "~+" and ingredient != "": + + operation = da.add if ingredient == "+" else da.subtract + instructions[recipe_index].append(operation) + continue + + component = model_components.get(ingredient) + + if component.type == "tigger-lsm": + + filename = component.path_or_name + tags = component.tags or () + + if not os.path.isfile(filename): + raise FileNotFoundError("{} not found.".format(filename)) + + sky_model = sky_model_nt(filename, tags) + sky_models.add(sky_model) + instructions[recipe_index].append(sky_model) + + elif component.type == "mds": + + filename = component.path_or_name + + options = ( + component.npix_x, + component.npix_y, + component.cellsize_x, + component.cellsize_y, + component.centre_x, + component.centre_y, + component.integrations_per_image, + component.channels_per_image, + ) + + if not os.path.exists(filename): + raise FileNotFoundError("{} not found.".format(filename)) + + degrid_model = degrid_model_nt(filename, *options) + degrid_models.add(degrid_model) + instructions[recipe_index].append(degrid_model) + + elif component.type == "column": + + column_name = component.path_or_name + + model_columns.add(column_name) + instructions[recipe_index].append(column_name) + + else: + instructions[recipe_index].append(ingredient) + + # TODO: Add message to log. + logger.info("The following model sources were obtained from " + "--input-model-recipe: \n" + " Columns: {} \n" + " Sky Models: {} \n" + " Degrid Models: {}", + model_columns or 'None', + {sm.name for sm in sky_models} or 'None', + {dm.name for dm in degrid_models} or 'None') + + # Generate a named tuple containing all the information required to + # build the model visibilities. + + model_recipe = Recipe( + Ingredients( + model_columns, + sky_models, + degrid_models + ), + instructions + ) + + if model_recipe.ingredients.sky_models: + logger.info("Recipe contains sky models - enabling prediction step.") + + if model_recipe.ingredients.degrid_models: + logger.info("Recipe contains degrid models - enabling degridding.") + + return model_recipe diff --git a/quartical/executor.py b/quartical/executor.py index ed0bcec7..1addfed7 100644 --- a/quartical/executor.py +++ b/quartical/executor.py @@ -44,6 +44,7 @@ def _execute(exitstack): output_opts = opts.output mad_flag_opts = opts.mad_flags dask_opts = opts.dask + model_components = internal.get_component_dict(opts) chain = internal.gains_to_chain(opts) # Special handling. # Init the logging proxy - an object which helps us ensure that logging @@ -60,7 +61,15 @@ def _execute(exitstack): # Now that we know where to put the log, log the final config state. parser.log_final_config(opts, config_files) - model_vis_recipe = preprocess.transcribe_recipe(model_opts.recipe) + # TODO: Deprecate legacy models. + if model_opts.recipe: + model_vis_recipe = preprocess.transcribe_legacy_recipe( + model_opts.recipe + ) + else: + model_vis_recipe = preprocess.transcribe_recipe( + model_opts.advanced_recipe, model_components + ) if dask_opts.scheduler == "distributed": From aed942a23a9aa0f81867f9dd103cd8aebf6372ff Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 15:17:44 +0200 Subject: [PATCH 09/63] Change advanced model to a simple boolean flag. --- quartical/config/argument_schema.yaml | 5 +++-- quartical/config/config_classes.py | 3 --- quartical/config/external.py | 7 ++++--- quartical/executor.py | 10 +++++----- 4 files changed, 12 insertions(+), 13 deletions(-) diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index 28644749..8deb2d72 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -120,9 +120,10 @@ input_model: (the default) will use an identity model. advanced_recipe: - dtype: Optional[str] + dtype: bool + default: false info: - Optional advanced recipe specification for use with degridder/more + Enable advanced recipe specification for use with degridder/more complicated modes of model construction. beam: diff --git a/quartical/config/config_classes.py b/quartical/config/config_classes.py index 9baffab9..d06b7d39 100644 --- a/quartical/config/config_classes.py +++ b/quartical/config/config_classes.py @@ -111,9 +111,6 @@ def __input_model_post_init__(self): self.__validate_choices__() self.__validate_element_choices__() - assert not (self.recipe and self.advanced_recipe), \ - "recipe and advanced_recipe are mutually exclusive." - def __model_component_post_init__(self): self.__validate_choices__() diff --git a/quartical/config/external.py b/quartical/config/external.py index 5c3de452..f3246d7b 100644 --- a/quartical/config/external.py +++ b/quartical/config/external.py @@ -24,13 +24,14 @@ def finalize_structure(additional_config): # Get last specified version of input_model.recipe. for cfg in additional_config[::-1]: - recipe = oc.select(cfg, "input_model.advanced_recipe") - if recipe is not None: + advanced_recipe = oc.select(cfg, "input_model.advanced_recipe") + recipe = oc.select(cfg, "input_model.recipe") + if recipe is not None and advanced_recipe: ingredients = re.split(r'([\+~:])', recipe) ingredients = [ i for i in ingredients if not bool(re.search(r'([\+~:])', i)) ] - models = set(i.split("@")[0] for i in ingredients) + models = list(dict.fromkeys(i.split("@")[0] for i in ingredients)) break FinalConfig = make_dataclass( diff --git a/quartical/executor.py b/quartical/executor.py index 1addfed7..69f91a11 100644 --- a/quartical/executor.py +++ b/quartical/executor.py @@ -62,13 +62,13 @@ def _execute(exitstack): parser.log_final_config(opts, config_files) # TODO: Deprecate legacy models. - if model_opts.recipe: - model_vis_recipe = preprocess.transcribe_legacy_recipe( - model_opts.recipe + if model_opts.advanced_recipe: + model_vis_recipe = preprocess.transcribe_recipe( + model_opts.recipe, model_components ) else: - model_vis_recipe = preprocess.transcribe_recipe( - model_opts.advanced_recipe, model_components + model_vis_recipe = preprocess.transcribe_legacy_recipe( + model_opts.recipe ) if dask_opts.scheduler == "distributed": From 0d43f3e386c2187777de59791f3fbccfb23e8109 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 16:16:06 +0200 Subject: [PATCH 10/63] Renaming. --- quartical/config/preprocess.py | 16 +++---- quartical/data_handling/degridder.py | 62 +++++++++++++++++----------- 2 files changed, 46 insertions(+), 32 deletions(-) diff --git a/quartical/config/preprocess.py b/quartical/config/preprocess.py index c016067d..125e86bb 100644 --- a/quartical/config/preprocess.py +++ b/quartical/config/preprocess.py @@ -14,14 +14,14 @@ "degrid_model_nt", ( "name", - "nxo", - "nyo", - "cellxo", - "cellyo", - "x0o", - "y0o", - "ipi", - "cpi" + "npix_x", + "npix_y", + "cellsize_x", + "cellsize_y", + "centre_x", + "centre_y", + "integrations_per_image", + "channels_per_image" ) ) diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py index 24711897..4852c3a6 100644 --- a/quartical/data_handling/degridder.py +++ b/quartical/data_handling/degridder.py @@ -10,14 +10,22 @@ from quartical.utils.collections import freeze_default_dict -def _degrid(time, freq, uvw, model): - - name, nxo, nyo, cellxo, cellyo, x0o, y0o, ipi, cpi = model +def _degrid(time, freq, uvw, component): + + name = component.name + npix_x = component.npix_x # Will be used in interpolation mode. + npix_y = component.npix_y # Will be used in interpolation mode. + cellsize_x = component.cellsize_x + cellsize_y = component.cellsize_y + centre_x = component.centre_x + centre_y = component.centre_y + integrations_per_image = component.integrations_per_image + channels_per_image = component.channels_per_image # TODO: Dodgy pattern, as this will produce dask arrays which will # need to be evaluated i.e. task-in-task which is a bad idea. OK for # the purposes of a prototype. - model_xds = xds_from_zarr(model.name)[0] + model_xds = xds_from_zarr(name)[0] # TODO: We want to go from the model to an image cube appropriate for this # chunks of data. The dimensions of the image cube will be determined @@ -26,29 +34,35 @@ def _degrid(time, freq, uvw, model): utime, utime_inv = np.unique(time, return_inverse=True) n_utime = utime.size - ipi = ipi or n_utime # Catch zero case. NOTE: -1? + ipi = integrations_per_image or n_utime n_mean_times = int(np.ceil(n_utime / ipi)) n_freq = freq.size + cpi = channels_per_image or n_freq n_mean_freqs = int(np.ceil(n_freq / cpi)) # Let's start off by simply reconstructing the model image as generated # by pfb-clean. - nxi, nyi = model_xds.npix_x, model_xds.npix_y - image_in = np.zeros((nxi, nyi), dtype=float) + native_npix_x = model_xds.npix_x + native_npix_y = model_xds.npix_y + native_image = np.zeros((native_npix_x, native_npix_y), dtype=float) + + # Sey up sympy symbols for expression evaluation. params = sm.symbols(('t', 'f')) params += sm.symbols(tuple(model_xds.params.values)) symexpr = parse_expr(model_xds.parametrisation) - modelf = lambdify(params, symexpr) - texpr = parse_expr(model_xds.texpr) - tfunc = lambdify(params[0], texpr) - fexpr = parse_expr(model_xds.fexpr) - ffunc = lambdify(params[1], fexpr) - Ix = model_xds.location_x.values - Iy = model_xds.location_y.values - coeffs = model_xds.coefficients.values + model_func = lambdify(params, symexpr) + time_expr = parse_expr(model_xds.texpr) + time_func = lambdify(params[0], time_expr) + freq_expr = parse_expr(model_xds.fexpr) + freq_func = lambdify(params[1], freq_expr) + + pixel_xs = model_xds.location_x.values + pixel_ys = model_xds.location_y.values + pixel_coeffs = model_xds.coefficients.values + # TODO: How do we handle the correlation axis neatly? vis = np.empty((time.size, freq.size, 4), dtype=np.complex128) for ti in range(n_mean_times): @@ -62,8 +76,8 @@ def _degrid(time, freq, uvw, model): degrid_time = utime[time_sel] mean_time = degrid_time.mean() - image_in[Ix, Iy] = modelf( - tfunc(mean_time), ffunc(mean_freq), *coeffs + native_image[pixel_xs, pixel_ys] = model_func( + time_func(mean_time), freq_func(mean_freq), *pixel_coeffs ) # NOTE: Select out appropriate rows for ipi and make selection @@ -75,19 +89,19 @@ def _degrid(time, freq, uvw, model): vis=vis[row_sel, freq_sel, 0], uvw=uvw, freq=degrid_freq, - dirty=image_in, - pixsize_x=model_xds.cell_rad_x, # Should be output value. - pixsize_y=model_xds.cell_rad_y, # Should be output value. - center_x=model_xds.center_x, # Should be output value. - center_y=model_xds.center_y, # Should be output value. + dirty=native_image, + pixsize_x=cellsize_x or model_xds.cell_rad_x, + pixsize_y=cellsize_y or model_xds.cell_rad_y, + center_x=centre_x or model_xds.center_x, + center_y=centre_y or model_xds.center_y, epsilon=1e-7, # TODO: Is this too high? do_wgridding=True, # Should be ok to leave True. divide_by_n=False, # Until otherwise informed. nthreads=6 # Should be equivalent to solver threads. ) - # Zero the image array between image slices. - image_in[:, :] = 0 + # Zero the image array between image slices as a precaution. + native_image[:, :] = 0 # Degridder only produces I - will need to be more sophisticated. vis[..., -1] = vis[..., 0] From 6aa67f28ee6d54ea5e5ea87caa8bfdfa83c37766 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 16:46:22 +0200 Subject: [PATCH 11/63] Fix tag format. --- quartical/config/model_component_schema.yaml | 8 ++++---- quartical/config/preprocess.py | 3 +-- quartical/data_handling/degridder.py | 4 ++-- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/quartical/config/model_component_schema.yaml b/quartical/config/model_component_schema.yaml index b1629c23..1ed5b6a6 100644 --- a/quartical/config/model_component_schema.yaml +++ b/quartical/config/model_component_schema.yaml @@ -17,10 +17,10 @@ model_component: dtype: Optional[List[str]] info: Tag for use in the predict. - region: - dtype: Optional[str] - info: - Name of region file to use in degridder. + # region: # Add when implemented. + # dtype: Optional[str] + # info: + # Name of region file to use in degridder. npix_x: dtype: Optional[int] info: diff --git a/quartical/config/preprocess.py b/quartical/config/preprocess.py index 125e86bb..f8960031 100644 --- a/quartical/config/preprocess.py +++ b/quartical/config/preprocess.py @@ -217,7 +217,7 @@ def transcribe_recipe(user_recipe, model_components): if not os.path.isfile(filename): raise FileNotFoundError("{} not found.".format(filename)) - sky_model = sky_model_nt(filename, tags) + sky_model = sky_model_nt(filename, tuple(tags)) sky_models.add(sky_model) instructions[recipe_index].append(sky_model) @@ -253,7 +253,6 @@ def transcribe_recipe(user_recipe, model_components): else: instructions[recipe_index].append(ingredient) - # TODO: Add message to log. logger.info("The following model sources were obtained from " "--input-model-recipe: \n" " Columns: {} \n" diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py index 4852c3a6..cdf0ee4b 100644 --- a/quartical/data_handling/degridder.py +++ b/quartical/data_handling/degridder.py @@ -97,13 +97,13 @@ def _degrid(time, freq, uvw, component): epsilon=1e-7, # TODO: Is this too high? do_wgridding=True, # Should be ok to leave True. divide_by_n=False, # Until otherwise informed. - nthreads=6 # Should be equivalent to solver threads. + nthreads=6 # TODO: Should be equivalent to solver threads. ) # Zero the image array between image slices as a precaution. native_image[:, :] = 0 - # Degridder only produces I - will need to be more sophisticated. + # TODO: Degridder only produces I - will need to be more sophisticated. vis[..., -1] = vis[..., 0] return vis From bba12a2fd91a8f44f6b3ddd9113f68f61a56b481 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 17:09:05 +0200 Subject: [PATCH 12/63] Make dask usage less unholy. --- quartical/data_handling/degridder.py | 44 +++++++++++++++------------- 1 file changed, 24 insertions(+), 20 deletions(-) diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py index cdf0ee4b..fefc869a 100644 --- a/quartical/data_handling/degridder.py +++ b/quartical/data_handling/degridder.py @@ -1,6 +1,7 @@ from collections import defaultdict import numpy as np import dask.array as da +import xarray from daskms.experimental.zarr import xds_from_zarr from scipy.interpolate import RegularGridInterpolator import sympy as sm @@ -10,9 +11,8 @@ from quartical.utils.collections import freeze_default_dict -def _degrid(time, freq, uvw, component): +def _degrid(time, freq, uvw, pixel_coeffs, component, meta_xds): - name = component.name npix_x = component.npix_x # Will be used in interpolation mode. npix_y = component.npix_y # Will be used in interpolation mode. cellsize_x = component.cellsize_x @@ -22,11 +22,6 @@ def _degrid(time, freq, uvw, component): integrations_per_image = component.integrations_per_image channels_per_image = component.channels_per_image - # TODO: Dodgy pattern, as this will produce dask arrays which will - # need to be evaluated i.e. task-in-task which is a bad idea. OK for - # the purposes of a prototype. - model_xds = xds_from_zarr(name)[0] - # TODO: We want to go from the model to an image cube appropriate for this # chunks of data. The dimensions of the image cube will be determined # by the time and freqency dimensions of the chunk in conjunction with @@ -44,23 +39,22 @@ def _degrid(time, freq, uvw, component): # Let's start off by simply reconstructing the model image as generated # by pfb-clean. - native_npix_x = model_xds.npix_x - native_npix_y = model_xds.npix_y + native_npix_x = meta_xds.npix_x + native_npix_y = meta_xds.npix_y native_image = np.zeros((native_npix_x, native_npix_y), dtype=float) # Sey up sympy symbols for expression evaluation. params = sm.symbols(('t', 'f')) - params += sm.symbols(tuple(model_xds.params.values)) - symexpr = parse_expr(model_xds.parametrisation) + params += sm.symbols(tuple(meta_xds.params.values)) + symexpr = parse_expr(meta_xds.parametrisation) model_func = lambdify(params, symexpr) - time_expr = parse_expr(model_xds.texpr) + time_expr = parse_expr(meta_xds.texpr) time_func = lambdify(params[0], time_expr) - freq_expr = parse_expr(model_xds.fexpr) + freq_expr = parse_expr(meta_xds.fexpr) freq_func = lambdify(params[1], freq_expr) - pixel_xs = model_xds.location_x.values - pixel_ys = model_xds.location_y.values - pixel_coeffs = model_xds.coefficients.values + pixel_xs = meta_xds.location_x.values + pixel_ys = meta_xds.location_y.values # TODO: How do we handle the correlation axis neatly? vis = np.empty((time.size, freq.size, 4), dtype=np.complex128) @@ -90,10 +84,10 @@ def _degrid(time, freq, uvw, component): uvw=uvw, freq=degrid_freq, dirty=native_image, - pixsize_x=cellsize_x or model_xds.cell_rad_x, - pixsize_y=cellsize_y or model_xds.cell_rad_y, - center_x=centre_x or model_xds.center_x, - center_y=centre_y or model_xds.center_y, + pixsize_x=cellsize_x or meta_xds.cell_rad_x, + pixsize_y=cellsize_y or meta_xds.cell_rad_y, + center_x=centre_x or meta_xds.center_x, + center_y=centre_y or meta_xds.center_y, epsilon=1e-7, # TODO: Is this too high? do_wgridding=True, # Should be ok to leave True. divide_by_n=False, # Until otherwise informed. @@ -191,12 +185,22 @@ def degrid(data_xds_list, model_vis_recipe, ms_path, model_opts): for degrid_model in degrid_models: + model_xds = xds_from_zarr(degrid_model.name)[0] + + # NOTE: This is convenient but will result in some extra + # information being embedded in the graph. Shouldn't be a problem. + meta_xds = xarray.Dataset( + coords=model_xds.coords, attrs=model_xds.attrs + ) + degrid_vis = da.blockwise( _degrid, ("rowlike", "chan", "corr"), data_xds.TIME.data, ("rowlike",), data_xds.CHAN_FREQ.data, ("chan",), data_xds.UVW.data, ("rowlike", "uvw"), + model_xds.coefficients.data, ("params", "comps"), degrid_model, None, + meta_xds, None, concatenate=True, align_arrays=False, meta=np.empty([0, 0, 0], dtype=np.complex128), From 808e17a35ce8e32b7b0dc496643d2e644b00da7c Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 17:28:55 +0200 Subject: [PATCH 13/63] Add clone. --- quartical/data_handling/degridder.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py index fefc869a..c325385e 100644 --- a/quartical/data_handling/degridder.py +++ b/quartical/data_handling/degridder.py @@ -1,6 +1,7 @@ from collections import defaultdict import numpy as np import dask.array as da +from dask.graph_manipulation import clone import xarray from daskms.experimental.zarr import xds_from_zarr from scipy.interpolate import RegularGridInterpolator @@ -196,7 +197,7 @@ def degrid(data_xds_list, model_vis_recipe, ms_path, model_opts): degrid_vis = da.blockwise( _degrid, ("rowlike", "chan", "corr"), data_xds.TIME.data, ("rowlike",), - data_xds.CHAN_FREQ.data, ("chan",), + clone(data_xds.CHAN_FREQ.data), ("chan",), data_xds.UVW.data, ("rowlike", "uvw"), model_xds.coefficients.data, ("params", "comps"), degrid_model, None, From 519106b88cfce0b21bf3b2207eaceb0081c05465 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 31 Jul 2023 17:58:23 +0200 Subject: [PATCH 14/63] Don't init vis as empty - duh. --- quartical/data_handling/degridder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py index c325385e..82c8463c 100644 --- a/quartical/data_handling/degridder.py +++ b/quartical/data_handling/degridder.py @@ -58,7 +58,7 @@ def _degrid(time, freq, uvw, pixel_coeffs, component, meta_xds): pixel_ys = meta_xds.location_y.values # TODO: How do we handle the correlation axis neatly? - vis = np.empty((time.size, freq.size, 4), dtype=np.complex128) + vis = np.zeros((time.size, freq.size, 4), dtype=np.complex128) for ti in range(n_mean_times): for fi in range(n_mean_freqs): From 58eff1927ca800e912e8434a29f64f0b10111c82 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Tue, 1 Aug 2023 11:40:10 +0200 Subject: [PATCH 15/63] Fix correlation handling in degrid code. --- quartical/data_handling/degridder.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py index 82c8463c..07aa2ed0 100644 --- a/quartical/data_handling/degridder.py +++ b/quartical/data_handling/degridder.py @@ -12,7 +12,7 @@ from quartical.utils.collections import freeze_default_dict -def _degrid(time, freq, uvw, pixel_coeffs, component, meta_xds): +def _degrid(time, freq, uvw, pixel_coeffs, component, meta_xds, n_corr): npix_x = component.npix_x # Will be used in interpolation mode. npix_y = component.npix_y # Will be used in interpolation mode. @@ -58,7 +58,7 @@ def _degrid(time, freq, uvw, pixel_coeffs, component, meta_xds): pixel_ys = meta_xds.location_y.values # TODO: How do we handle the correlation axis neatly? - vis = np.zeros((time.size, freq.size, 4), dtype=np.complex128) + vis = np.zeros((time.size, freq.size, n_corr), dtype=np.complex128) for ti in range(n_mean_times): for fi in range(n_mean_freqs): @@ -182,6 +182,8 @@ def degrid(data_xds_list, model_vis_recipe, ms_path, model_opts): for data_xds in data_xds_list: + n_corr = data_xds.dims["corr"] + model_vis = defaultdict(list) for degrid_model in degrid_models: @@ -202,10 +204,11 @@ def degrid(data_xds_list, model_vis_recipe, ms_path, model_opts): model_xds.coefficients.data, ("params", "comps"), degrid_model, None, meta_xds, None, + n_corr, None, concatenate=True, align_arrays=False, meta=np.empty([0, 0, 0], dtype=np.complex128), - new_axes={"corr": 4} # TODO: Shouldn't be hardcoded. + new_axes={"corr": n_corr} ) model_vis[degrid_model].append(degrid_vis) From 2371715cfb124a867e58c016eda6e8ea954103ec Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Tue, 1 Aug 2023 11:54:12 +0200 Subject: [PATCH 16/63] Expose degrid (predict) threading in model arguments. --- quartical/config/argument_schema.yaml | 9 +++++++++ quartical/data_handling/degridder.py | 16 ++++++++++++++-- 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index 8deb2d72..ac7d768c 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -126,6 +126,15 @@ input_model: Enable advanced recipe specification for use with degridder/more complicated modes of model construction. + threads: + dtype: int + default: 1 + info: + Controls the number of threads used internally by the various predict + methods. This should typically be set to the same value as + solver.threads if solver.threads is in use. Currently only supported + for degridding. + beam: dtype: Optional[str] info: diff --git a/quartical/data_handling/degridder.py b/quartical/data_handling/degridder.py index 07aa2ed0..f4630b37 100644 --- a/quartical/data_handling/degridder.py +++ b/quartical/data_handling/degridder.py @@ -12,7 +12,16 @@ from quartical.utils.collections import freeze_default_dict -def _degrid(time, freq, uvw, pixel_coeffs, component, meta_xds, n_corr): +def _degrid( + time, + freq, + uvw, + pixel_coeffs, + component, + meta_xds, + n_corr, + n_thread +): npix_x = component.npix_x # Will be used in interpolation mode. npix_y = component.npix_y # Will be used in interpolation mode. @@ -92,7 +101,7 @@ def _degrid(time, freq, uvw, pixel_coeffs, component, meta_xds, n_corr): epsilon=1e-7, # TODO: Is this too high? do_wgridding=True, # Should be ok to leave True. divide_by_n=False, # Until otherwise informed. - nthreads=6 # TODO: Should be equivalent to solver threads. + nthreads=n_thread ) # Zero the image array between image slices as a precaution. @@ -178,6 +187,8 @@ def degrid(data_xds_list, model_vis_recipe, ms_path, model_opts): degrid_models = model_vis_recipe.ingredients.degrid_models + n_thread = model_opts.threads + degrid_list = [] for data_xds in data_xds_list: @@ -205,6 +216,7 @@ def degrid(data_xds_list, model_vis_recipe, ms_path, model_opts): degrid_model, None, meta_xds, None, n_corr, None, + n_thread, None, concatenate=True, align_arrays=False, meta=np.empty([0, 0, 0], dtype=np.complex128), From 7f1bc9aa7dab125474275c125bdc3f7d01465dce Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Tue, 1 Aug 2023 12:42:11 +0200 Subject: [PATCH 17/63] Disable degrid inputs in non-advanced mode. --- quartical/config/preprocess.py | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/quartical/config/preprocess.py b/quartical/config/preprocess.py index f8960031..13dd69be 100644 --- a/quartical/config/preprocess.py +++ b/quartical/config/preprocess.py @@ -6,7 +6,6 @@ import os.path from dataclasses import dataclass from typing import List, Dict, Set, Any -from ast import literal_eval sky_model_nt = namedtuple("sky_model_nt", ("name", "tags")) @@ -108,15 +107,10 @@ def transcribe_legacy_recipe(user_recipe): elif ".mds" in ingredient: - filename, _, options = ingredient.partition("@") - options = literal_eval(options) # Add fail on missing option. - - if not os.path.exists(filename): - raise FileNotFoundError("{} not found.".format(filename)) - - degrid_model = degrid_model_nt(filename, *options) - degrid_models.add(degrid_model) - instructions[recipe_index].append(degrid_model) + # TODO: Add link to documentation. + raise ValueError( + ".mds inputs are only supported in advanced model mode." + ) elif ingredient != "": model_columns.add(ingredient) @@ -150,9 +144,6 @@ def transcribe_legacy_recipe(user_recipe): if model_recipe.ingredients.sky_models: logger.info("Recipe contains sky models - enabling prediction step.") - if model_recipe.ingredients.degrid_models: - logger.info("Recipe contains degrid models - enabling degridding.") - return model_recipe From 00af5188547223cf8391f47b5eaeac98fc4c821f Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Tue, 1 Aug 2023 15:45:34 +0200 Subject: [PATCH 18/63] Use legacy models in tests until such time as I can reconfigure them. --- testing/fixtures/config.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testing/fixtures/config.py b/testing/fixtures/config.py index 7850157c..bbc9830a 100644 --- a/testing/fixtures/config.py +++ b/testing/fixtures/config.py @@ -1,5 +1,5 @@ import pytest -from quartical.config.preprocess import transcribe_recipe +from quartical.config.preprocess import transcribe_legacy_recipe from quartical.config.internal import gains_to_chain @@ -40,4 +40,4 @@ def chain(opts): @pytest.fixture(scope="module") def recipe(model_opts): - return transcribe_recipe(model_opts.recipe) + return transcribe_legacy_recipe(model_opts.recipe) From e3340264773d53860a7b17eaf793d063de73cbdc Mon Sep 17 00:00:00 2001 From: landmanbester Date: Wed, 2 Aug 2023 15:50:50 +0200 Subject: [PATCH 19/63] inconsistent RNG error --- quartical/config/gain_schema.yaml | 2 +- quartical/gains/complex/__init__.py | 4 + quartical/interpolation/interpolants.py | 124 ++++++++++++++++++++++++ quartical/interpolation/interpolate.py | 93 ++++++++++++++---- 4 files changed, 203 insertions(+), 20 deletions(-) diff --git a/quartical/config/gain_schema.yaml b/quartical/config/gain_schema.yaml index 0a462bb8..070028e2 100644 --- a/quartical/config/gain_schema.yaml +++ b/quartical/config/gain_schema.yaml @@ -71,7 +71,7 @@ gain: interp_method: dtype: str default: "2dlinear" - choices: ["2dlinear", "2dspline"] + choices: ["2dlinear", "2dspline", "1dsmooth"] info: Set interpolation method. diff --git a/quartical/gains/complex/__init__.py b/quartical/gains/complex/__init__.py index 39e6cd96..d0c40e89 100644 --- a/quartical/gains/complex/__init__.py +++ b/quartical/gains/complex/__init__.py @@ -34,3 +34,7 @@ class DiagComplex(Complex): def __init__(self, term_name, term_opts): super().__init__(term_name, term_opts) + + # needed for smoothing + if 'jhj' not in self.interpolation_targets: + self.interpolation_targets.append('jhj') diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 0c384645..5151ed77 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -402,3 +402,127 @@ def interpolate_missing(input_xda, mode="normal"): ) return input_xda.copy(data=interp) + + +def smooth_ampphase(gains, + jhj, + flags, + ti, + fi, + to, + fo, + combine_by_time=True, + niter=5, + nu0=2.0, + padding=1.2): + import nifty8 as ift + from ducc0.fft import good_size + if not combine_by_time: + raise NotImplementedError + # weighted sum over time axis + ntimei, nchani, nanti, ndiri, ncorri = gains.shape + ntimeo = to.size + nfreqo = fo.size + assert nanti == 1 + assert ndiri == 1 + assert ncorri == 1 + gains = gains.squeeze() + jhj = jhj.squeeze().real + gain = np.zeros((nchani), dtype=gains.dtype) + wgt = np.zeros((nchani), dtype=jhj.dtype) + flags = flags.squeeze().astype(bool) + jhj[flags] = 0.0 + for t in range(ntimei): + gain += jhj[t].squeeze() * gains[t].squeeze() + wgt += jhj[t].squeeze() + # normalise by sum of weights + mask = wgt > 0 + if not mask.any(): + return np.ones((ntimeo, nfreqo, nanti, ndiri, ncorri), + dtype=gains.dtype) + + print(wgt[mask].max(), wgt[mask].min()) + gain[mask] = gain[mask]/wgt[mask] + + # set up correlated field model (hardcoding hypers for now) + npix_f = good_size(int(nchani*padding)) + pospace_large = ift.RGSpace([npix_f]) + pospace_small = ift.RGSpace([nchani]) + spfreq_amp = ift.RGSpace(npix_f) + spfreq_phase = ift.RGSpace(npix_f) + + # amplitude model + cfmakera = ift.CorrelatedFieldMaker('amplitude') + cfmakera.add_fluctuations(spfreq_amp, (0.1, 1e-2), None, None, (-3, 1), + 'f') + cfmakera.set_amplitude_total_offset(0., (1e-2, 1e-6)) + cf_amp = cfmakera.finalize() + + normalized_amp = cfmakera.get_normalized_amplitudes() + pspec_freq_amp = normalized_amp[0]**2 + + # phase model + cfmakerp = ift.CorrelatedFieldMaker('phase') + cfmakerp.add_fluctuations(spfreq_phase, (0.1, 1e-2), None, None, (-3, 1), + 'f') + cfmakerp.set_amplitude_total_offset(0., (1e-2, 1e-6)) + cf_phase = cfmakerp.finalize() + + normalized_phase = cfmakerp.get_normalized_amplitudes() + pspec_freq_phase = normalized_phase[0]**2 + + TF = ift.SliceOperator(cf_amp.target, pospace_small.shape) + + signal = TF @ ift.exp(cf_amp.real + 1j*cf_phase.real) + + cf_amp_small = TF @ cf_amp + cf_phase_small = TF @ cf_phase + + tmp = ift.makeField(signal.target, ~mask) + # TODO - this should include averaging + R = ift.MaskOperator(tmp) + dspace = R.target + data = ift.makeField(dspace, gain[mask]) + signal_response = R(signal) + + # Minimization parameters + n_samples = 3 + n_iterations = 5 + ic_sampling = ift.AbsDeltaEnergyController(name='Sampling', deltaE=0.01, iteration_limit=50) + ic_newton = ift.AbsDeltaEnergyController(name='Newton', deltaE=0.01, iteration_limit=15) + minimizer = ift.NewtonCG(ic_newton, enable_logging=True) + + # Set up likelihood energy and information Hamiltonian + W = wgt[mask] + assert (W > 0).all() + covinv = ift.makeField(dspace, W) + Ninv = ift.DiagonalOperator(covinv, sampling_dtype=complex) + inp = (ift.Adder(data, neg=True) @ signal_response) # required to compute the resid for VariableCovGaussianEnergy + BC = ift.ContractionOperator(dspace, None).adjoint # broadcast scalar over everything (None) + weightop = ift.DiagonalOperator(ift.makeField(dspace, W)) + scalarop = ift.GammaOperator(BC.domain, mean=2.0, var=2).ducktape('icov_scalar') + icovop = weightop @ BC @ scalarop.real + inp = inp.ducktape_left('resid') + icovop.ducktape_left('icov') + likelihood_energy = ift.VariableCovarianceGaussianEnergy(R.target, 'resid', 'icov', + data.dtype) @ inp + + samples = ift.optimize_kl(likelihood_energy, + n_iterations, + n_samples, + minimizer, + ic_sampling, + None) # for GeoVI + + # get the mean + signal_mean = samples.average(signal.force).val + + # interpolate to output frequencies + amp = np.interp(fo, fi, np.abs(signal_mean)) + phase = np.interp(fo, fi, np.angle(signal_mean)) + signal_mean = amp*np.exp(1j*phase) + + # broadcast to output times + signal_mean = np.tile(signal_mean[None, :], (to.size, 1)) + + # broadcast to expected output shape + return signal_mean[:, :, None, None, None] diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index b37b39e1..d7be4ae9 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -92,25 +92,33 @@ def load_and_interpolate_gains(gain_xds_lod, chain, output_directory): # Remove time/chan chunking and rechunk by antenna. merged_xds = merged_xds.chunk({**merged_xds.dims, "antenna": 1}) - # Create a converter object to handle moving between native and - # interpolation representations. - converter = Converter(term) - - # Convert standard representation to a representation which can - # be interpolated. Replace flags with NaNs. - merged_xds = convert_native_to_interp(merged_xds, converter) - - # Interpolate onto the given grids. TODO: Add back spline support. - interpolated_xds_list = [ - interpolate(merged_xds, xds, term) for xds in term_xds_list - ] - - # Convert from representation which can be interpolated back into - # native representation. - interpolated_xds_list = [ - convert_interp_to_native(xds, converter) - for xds in interpolated_xds_list - ] + if term.interp_method == "1dsmooth": + if merged_xds.TYPE != "diag_complex": + raise ValueError("Smoothing only supported for diag_complex type") + ntime, nchan, nant, ndir, ncorr = merged_xds.gains.shape + if ndir > 1: + raise ValueError("Smoothing only supported for direction independent gains") + interpolated_xds_list = bsmooth(merged_xds, term_xds_list) + else: + # Create a converter object to handle moving between native and + # interpolation representations. + converter = Converter(term) + + # Convert standard representation to a representation which can + # be interpolated. Replace flags with NaNs. + merged_xds = convert_native_to_interp(merged_xds, converter) + + # Interpolate onto the given grids. TODO: Add back spline support. + interpolated_xds_list = [ + interpolate(merged_xds, xds, term) for xds in term_xds_list + ] + + # Convert from representation which can be interpolated back into + # native representation. + interpolated_xds_list = [ + convert_interp_to_native(xds, converter) + for xds in interpolated_xds_list + ] # Make the interpolated xds consistent with the current run. interpolated_xds_list = [ @@ -275,3 +283,50 @@ def compute_and_reload(directory, gain_xds_list): rxds[k] = lxds[k] return gain_xds_list + + +def bsmooth(merged_xds, target_xds): + out_time = [] + out_freq =[] + for ds in target_xds: + out_time.append(ds.gain_time.data) + out_freq.append(ds.gain_freq.data) + out_time = np.unique(np.concatenate(out_time)) + out_freq = np.unique(np.concatenate(out_freq)) + ntimeo = out_time.size + nfreqo = out_freq.size + ntimei, nfreqi, nant, ndir, ncorr = merged_xds.gains.shape + + # we also want to chunk by correlation since they can be smoothed separately + merged_xds = merged_xds.chunk({"correlation": 1}) + + from quartical.interpolation.interpolants import smooth_ampphase + combine_by_time = True + smoothed_gains = da.blockwise( + smooth_ampphase, 'tfpdc', + merged_xds.gains.data, 'tfpdc', + merged_xds.jhj.data, 'tfpdc', + merged_xds.gain_flags.data, 'tfpd', + merged_xds.gain_time.data, None, + merged_xds.gain_freq.data, None, + out_time, None, + out_freq, None, + combine_by_time, None, + dtype=merged_xds.gains.dtype, + align_arrays=False, + adjust_chunks={'t': ntimeo, 'f': nfreqo}, + meta=np.empty((ntimeo, nfreqo, nant, ndir, ncorr), + dtype=merged_xds.gains.dtype) + ) + + for i, ds in enumerate(target_xds): + t = ds.gain_time.data + _, idx0, idx1 = np.intersect1d(t, out_time, + assume_unique=True, + return_indices=True) + ds['gains'] = (ds.GAIN_AXES, smoothed_gains[idx1]) + + target_xds = [ds.chunk({"correlation": ncorr}) for ds in target_xds] + + return target_xds + From aaaf76afb3978c2f292b1ab6d36263d648593e40 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Thu, 3 Aug 2023 18:35:25 +0200 Subject: [PATCH 20/63] Utilise environment variable when dask.address is unset. (#288) * Fix version drift. * Bump to 0.2.0 * Inspect envvar for scheduler address when one isn't specified. * Encode environment varraible as ascii. * Simplify. --- quartical/executor.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/quartical/executor.py b/quartical/executor.py index ed0bcec7..06fc85ae 100644 --- a/quartical/executor.py +++ b/quartical/executor.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- from contextlib import ExitStack +import os import sys from loguru import logger import dask @@ -68,9 +69,14 @@ def _execute(exitstack): # distributed enviroment. This *may* be dangerous. Monitor. dask.config.set({"distributed.worker.daemon": False}) - if dask_opts.address: - logger.info("Initializing distributed client.") - client = exitstack.enter_context(Client(dask_opts.address)) + address = dask_opts.address or os.environ.get("DASK_SCHEDULER_ADDRESS") + + if address: + logger.info( + f"Initializing distributed client using scheduler address: " + f"{address}" + ) + client = exitstack.enter_context(Client(address)) else: logger.info("Initializing distributed client using LocalCluster.") cluster = LocalCluster( From ad4f91bb2869eeedb975a8d5baed3e4ebdbbf5b0 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 4 Aug 2023 11:31:23 +0200 Subject: [PATCH 21/63] Fix incorrect import in tests. --- testing/tests/config/test_preprocess.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/testing/tests/config/test_preprocess.py b/testing/tests/config/test_preprocess.py index ed6032d3..6b983725 100644 --- a/testing/tests/config/test_preprocess.py +++ b/testing/tests/config/test_preprocess.py @@ -1,5 +1,5 @@ import pytest -from quartical.config.preprocess import transcribe_recipe, sky_model_nt +from quartical.config.preprocess import transcribe_legacy_recipe, sky_model_nt import dask.array as da import os.path @@ -51,7 +51,7 @@ # we do not attempt to validate column names in the preprocess step. invalid_recipes = { - "": ValueError, + # "": ValueError, # NOTE: This case may not be needed. Omitting. "dummy.lsm.html": FileNotFoundError } @@ -73,7 +73,7 @@ def test_transcribe_recipe_valid(valid_recipe, monkeypatch): # Patch isfile functionality to allow use of ficticious files. monkeypatch.setattr(os.path, "isfile", lambda filename: True) - recipe = transcribe_recipe(input_recipe) + recipe = transcribe_legacy_recipe(input_recipe) # Check that the opts has been updated with the correct internal recipe. assert recipe.instructions == expected_output @@ -87,4 +87,4 @@ def test_transcribe_recipe_invalid(invalid_recipe): input_recipe, expected_output = invalid_recipe with pytest.raises(expected_output): - transcribe_recipe(input_recipe) + transcribe_legacy_recipe(input_recipe) From 232ae41d937491b3d88a0134a6bf5f3595fa67e2 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 4 Aug 2023 16:00:55 +0200 Subject: [PATCH 22/63] Commit WIP code for experimental fragment branch. --- pyproject.toml | 2 +- quartical/apps/summary.py | 31 +++++++++++++++------------ quartical/data_handling/angles.py | 8 +++---- quartical/data_handling/chunking.py | 9 +++++--- quartical/data_handling/ms_handler.py | 29 +++++++++++++++++-------- quartical/data_handling/predict.py | 20 ++++++++++------- 6 files changed, 60 insertions(+), 39 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ceee48d2..df3ca683 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ columnar = "^1.4.1" "ruamel.yaml" = "^0.17.26" dask = {extras = ["diagnostics"], version = "^2023.1.0"} distributed = "^2023.1.0" -dask-ms = {extras = ["s3", "xarray", "zarr"], version = "^0.2.16"} +dask-ms = {git = "https://github.com/ratt-ru/dask-ms.git", branch = "multisource-experimental", extras = ["s3", "xarray", "zarr"]} codex-africanus = {extras = ["dask", "scipy", "astropy", "python-casacore"], version = "^0.3.4"} astro-tigger-lsm = "^1.7.2" loguru = "^0.7.0" diff --git a/quartical/apps/summary.py b/quartical/apps/summary.py index c6e5d915..1d0a8dd2 100644 --- a/quartical/apps/summary.py +++ b/quartical/apps/summary.py @@ -1,6 +1,9 @@ import argparse from pathlib import Path -from daskms import xds_from_storage_ms, xds_from_storage_table +from daskms.experimental.multisource import ( + xds_from_ms_fragment, + xds_from_table_fragment +) from daskms.fsspec_store import DaskMSStore import numpy as np import dask.array as da @@ -44,7 +47,7 @@ def configure_loguru(output_dir): def antenna_info(path): # NOTE: Assume one dataset for now. - ant_xds = xds_from_storage_table(path + "::ANTENNA")[0] + ant_xds = xds_from_table_fragment(path + "::ANTENNA")[0] antenna_names = ant_xds.NAME.values antenna_mounts = ant_xds.MOUNT.values @@ -64,7 +67,7 @@ def antenna_info(path): def data_desc_info(path): - dd_xds_list = xds_from_storage_table( # noqa + dd_xds_list = xds_from_table_fragment( # noqa path + "::DATA_DESCRIPTION", group_cols=["__row__"], chunks={"row": 1, "chan": -1} @@ -76,7 +79,7 @@ def data_desc_info(path): def feed_info(path): - feed_xds_list = xds_from_storage_table( + feed_xds_list = xds_from_table_fragment( path + "::FEED", group_cols=["SPECTRAL_WINDOW_ID"], chunks={"row": -1} @@ -106,7 +109,7 @@ def feed_info(path): def flag_cmd_info(path): - flag_cmd_xds = xds_from_storage_table(path + "::FLAG_CMD") # noqa + flag_cmd_xds = xds_from_table_fragment(path + "::FLAG_CMD") # noqa # Not printing any summary information for this subtable yet - not sure # what is relevant. @@ -114,7 +117,7 @@ def flag_cmd_info(path): def field_info(path): - field_xds = xds_from_storage_table(path + "::FIELD")[0] + field_xds = xds_from_table_fragment(path + "::FIELD")[0] ids = [i for i in field_xds.SOURCE_ID.values] names = [n for n in field_xds.NAME.values] @@ -141,7 +144,7 @@ def field_info(path): def history_info(path): - history_xds = xds_from_storage_table(path + "::HISTORY")[0] # noqa + history_xds = xds_from_table_fragment(path + "::HISTORY")[0] # noqa # Not printing any summary information for this subtable yet - not sure # what is relevant. @@ -149,7 +152,7 @@ def history_info(path): def observation_info(path): - observation_xds = xds_from_storage_table(path + "::OBSERVATION")[0] # noqa + observation_xds = xds_from_table_fragment(path + "::OBSERVATION")[0] # noqa # Not printing any summary information for this subtable yet - not sure # what is relevant. @@ -157,7 +160,7 @@ def observation_info(path): def polarization_info(path): - polarization_xds = xds_from_storage_table(path + "::POLARIZATION")[0] + polarization_xds = xds_from_table_fragment(path + "::POLARIZATION")[0] corr_types = polarization_xds.CORR_TYPE.values @@ -175,7 +178,7 @@ def polarization_info(path): def processor_info(path): - processor_xds = xds_from_storage_table(path + "::PROCESSOR")[0] # noqa + processor_xds = xds_from_table_fragment(path + "::PROCESSOR")[0] # noqa # Not printing any summary information for this subtable yet - not sure # what is relevant. @@ -183,7 +186,7 @@ def processor_info(path): def spw_info(path): - spw_xds_list = xds_from_storage_table( + spw_xds_list = xds_from_table_fragment( path + "::SPECTRAL_WINDOW", group_cols=["__row__"], chunks={"row": 1, "chan": -1} @@ -207,7 +210,7 @@ def spw_info(path): def state_info(path): - state_xds = xds_from_storage_table(path + "::STATE")[0] # noqa + state_xds = xds_from_table_fragment(path + "::STATE")[0] # noqa # Not printing any summary information for this subtable yet - not sure # what is relevant. @@ -226,7 +229,7 @@ def source_info(path): def pointing_info(path): - pointing_xds = xds_from_storage_table(path + "::POINTING")[0] # noqa + pointing_xds = xds_from_table_fragment(path + "::POINTING")[0] # noqa # Not printing any summary information for this subtable yet - not sure # what is relevant. @@ -355,7 +358,7 @@ def summary(): # Open the data, grouping by the usual columns. Use these datasets to # produce some useful summaries. - data_xds_list = xds_from_storage_ms( + data_xds_list = xds_from_ms_fragment( path, index_cols=("TIME",), columns=("TIME", "FLAG", "FLAG_ROW", "DATA"), diff --git a/quartical/data_handling/angles.py b/quartical/data_handling/angles.py index 343ab466..f7e6098e 100644 --- a/quartical/data_handling/angles.py +++ b/quartical/data_handling/angles.py @@ -2,7 +2,7 @@ import casacore.measures import casacore.quanta as pq -from daskms import xds_from_storage_table +from daskms.experimental.multisource import xds_from_table_fragment import dask.array as da import threading from dask.graph_manipulation import clone @@ -24,9 +24,9 @@ def make_parangle_xds_list(ms_path, data_xds_list): # This may need to be more sophisticated. TODO: Can we guarantee that # these only ever have one element? - anttab = xds_from_storage_table(ms_path + "::ANTENNA")[0] - feedtab = xds_from_storage_table(ms_path + "::FEED")[0] - fieldtab = xds_from_storage_table(ms_path + "::FIELD")[0] + anttab = xds_from_table_fragment(ms_path + "::ANTENNA")[0] + feedtab = xds_from_table_fragment(ms_path + "::FEED")[0] + fieldtab = xds_from_table_fragment(ms_path + "::FIELD")[0] # We do this eagerly to make life easier. feeds = feedtab.POLARIZATION_TYPE.values diff --git a/quartical/data_handling/chunking.py b/quartical/data_handling/chunking.py index ce6fe353..4003e737 100644 --- a/quartical/data_handling/chunking.py +++ b/quartical/data_handling/chunking.py @@ -1,7 +1,10 @@ import dask.delayed as dd import numpy as np import dask.array as da -from daskms import xds_from_storage_ms, xds_from_storage_table +from daskms.experimental.multisource import ( + xds_from_ms_fragment, + xds_from_table_fragment +) def compute_chunking(ms_opts, compute=True): @@ -10,7 +13,7 @@ def compute_chunking(ms_opts, compute=True): # necessary to determine initial chunking over row and chan. TODO: Test # multi-SPW/field cases. Implement a memory budget. - indexing_xds_list = xds_from_storage_ms( + indexing_xds_list = xds_from_ms_fragment( ms_opts.path, columns=("TIME", "INTERVAL"), index_cols=("TIME",), @@ -24,7 +27,7 @@ def compute_chunking(ms_opts, compute=True): compute=False ) - spw_xds_list = xds_from_storage_table( + spw_xds_list = xds_from_table_fragment( ms_opts.path + "::SPECTRAL_WINDOW", group_cols=["__row__"], columns=["CHAN_FREQ", "CHAN_WIDTH"], diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index 3b923432..edafa2fa 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -2,9 +2,12 @@ import warnings import dask.array as da import numpy as np -from daskms import (xds_from_storage_ms, - xds_from_storage_table, - xds_to_storage_table) +from daskms import xds_to_storage_table +from daskms.experimental.multisource import ( + xds_from_ms_fragment, + xds_from_table_fragment, + xds_to_table_fragment +) from dask.graph_manipulation import clone from loguru import logger from quartical.weights.weights import initialize_weights @@ -28,7 +31,7 @@ def read_xds_list(model_columns, ms_opts): data_xds_list: A list of appropriately chunked xarray datasets. """ - antenna_xds = xds_from_storage_table(ms_opts.path + "::ANTENNA")[0] + antenna_xds = xds_from_table_fragment(ms_opts.path + "::ANTENNA")[0] n_ant = antenna_xds.dims["row"] @@ -36,7 +39,7 @@ def read_xds_list(model_columns, ms_opts): "observation.", n_ant) # Determine the number/type of correlations present in the measurement set. - pol_xds = xds_from_storage_table(ms_opts.path + "::POLARIZATION")[0] + pol_xds = xds_from_table_fragment(ms_opts.path + "::POLARIZATION")[0] try: corr_types = [CORR_TYPES[ct] for ct in pol_xds.CORR_TYPE.values[0]] @@ -56,7 +59,7 @@ def read_xds_list(model_columns, ms_opts): # probably need to be done on a per xds basis. Can probably be accomplished # by merging the field xds grouped by DDID into data grouped by DDID. - field_xds = xds_from_storage_table(ms_opts.path + "::FIELD")[0] + field_xds = xds_from_table_fragment(ms_opts.path + "::FIELD")[0] phase_dir = np.squeeze(field_xds.PHASE_DIR.values) field_names = field_xds.NAME.values @@ -90,7 +93,7 @@ def read_xds_list(model_columns, ms_opts): schema[ms_opts.weight_column] = {'dims': ('chan', 'corr')} try: - data_xds_list = xds_from_storage_ms( + data_xds_list = xds_from_ms_fragment( ms_opts.path, columns=columns, index_cols=("TIME",), @@ -103,7 +106,7 @@ def read_xds_list(model_columns, ms_opts): f"Invalid/missing column specified. Underlying error: {e}." ) from e - spw_xds_list = xds_from_storage_table( + spw_xds_list = xds_from_table_fragment( ms_opts.path + "::SPECTRAL_WINDOW", group_cols=["__row__"], columns=["CHAN_FREQ", "CHAN_WIDTH"], @@ -213,7 +216,7 @@ def write_xds_list(xds_list, ref_xds_list, ms_path, output_opts): if not (output_opts.products or output_opts.flags): return [None] * len(xds_list) # Write nothing to the MS. - pol_xds = xds_from_storage_table(ms_path + "::POLARIZATION")[0] + pol_xds = xds_from_table_fragment(ms_path + "::POLARIZATION")[0] corr_types = [CORR_TYPES[ct] for ct in pol_xds.CORR_TYPE.values[0]] ms_n_corr = len(corr_types) @@ -302,6 +305,14 @@ def write_xds_list(xds_list, ref_xds_list, ms_path, output_opts): rechunk=True # Needed to ensure zarr chunks map correctly to disk. ) + # write_xds_list = xds_to_table_fragment( + # xds_list, + # "delta1.ms", + # ms_path, + # columns=output_cols, + # rechunk=True # Needed to ensure zarr chunks map correctly to disk. + # ) + return write_xds_list diff --git a/quartical/data_handling/predict.py b/quartical/data_handling/predict.py index 6bac1ef4..cc8bda04 100644 --- a/quartical/data_handling/predict.py +++ b/quartical/data_handling/predict.py @@ -7,7 +7,8 @@ import dask from xarray import DataArray, Dataset from dask.graph_manipulation import clone -from daskms import xds_from_storage_table +from daskms.experimental.multisource import xds_from_table_fragment + from loguru import logger import numpy as np import Tigger @@ -310,21 +311,24 @@ def get_support_tables(ms_path): "SPECTRAL_WINDOW", "POLARIZATION", "FEED")} # All rows at once - lazy_tables = {"ANTENNA": xds_from_storage_table(n["ANTENNA"]), - "FEED": xds_from_storage_table(n["FEED"])} + lazy_tables = {"ANTENNA": xds_from_table_fragment(n["ANTENNA"]), + "FEED": xds_from_table_fragment(n["FEED"])} compute_tables = { # NOTE: Even though this has a fixed shape, I have ammended it to # also group by row. This just makes life fractionally easier. - "DATA_DESCRIPTION": xds_from_storage_table(n["DATA_DESCRIPTION"], - group_cols="__row__"), + "DATA_DESCRIPTION": xds_from_table_fragment( + n["DATA_DESCRIPTION"], group_cols="__row__" + ), # Variably shaped, need a dataset per row "FIELD": - xds_from_storage_table(n["FIELD"], group_cols="__row__"), + xds_from_table_fragment(n["FIELD"], group_cols="__row__"), "SPECTRAL_WINDOW": - xds_from_storage_table(n["SPECTRAL_WINDOW"], group_cols="__row__"), + xds_from_table_fragment( + n["SPECTRAL_WINDOW"], group_cols="__row__" + ), "POLARIZATION": - xds_from_storage_table(n["POLARIZATION"], group_cols="__row__"), + xds_from_table_fragment(n["POLARIZATION"], group_cols="__row__"), } lazy_tables.update(dask.compute(compute_tables)[0]) From ec2f4967ea8e3ed5d0814b0e9b2db075d9dffe50 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 4 Aug 2023 16:24:14 +0200 Subject: [PATCH 23/63] Import from fragments. --- quartical/apps/summary.py | 2 +- quartical/data_handling/angles.py | 2 +- quartical/data_handling/chunking.py | 2 +- quartical/data_handling/ms_handler.py | 3 ++- quartical/data_handling/predict.py | 2 +- 5 files changed, 6 insertions(+), 5 deletions(-) diff --git a/quartical/apps/summary.py b/quartical/apps/summary.py index 1d0a8dd2..0d3dab90 100644 --- a/quartical/apps/summary.py +++ b/quartical/apps/summary.py @@ -1,6 +1,6 @@ import argparse from pathlib import Path -from daskms.experimental.multisource import ( +from daskms.experimental.fragments import ( xds_from_ms_fragment, xds_from_table_fragment ) diff --git a/quartical/data_handling/angles.py b/quartical/data_handling/angles.py index f7e6098e..7911a24f 100644 --- a/quartical/data_handling/angles.py +++ b/quartical/data_handling/angles.py @@ -2,7 +2,7 @@ import casacore.measures import casacore.quanta as pq -from daskms.experimental.multisource import xds_from_table_fragment +from daskms.experimental.fragments import xds_from_table_fragment import dask.array as da import threading from dask.graph_manipulation import clone diff --git a/quartical/data_handling/chunking.py b/quartical/data_handling/chunking.py index 4003e737..867a28a7 100644 --- a/quartical/data_handling/chunking.py +++ b/quartical/data_handling/chunking.py @@ -1,7 +1,7 @@ import dask.delayed as dd import numpy as np import dask.array as da -from daskms.experimental.multisource import ( +from daskms.experimental.fragments import ( xds_from_ms_fragment, xds_from_table_fragment ) diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index edafa2fa..e1682cfd 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -3,7 +3,7 @@ import dask.array as da import numpy as np from daskms import xds_to_storage_table -from daskms.experimental.multisource import ( +from daskms.experimental.fragments import ( xds_from_ms_fragment, xds_from_table_fragment, xds_to_table_fragment @@ -305,6 +305,7 @@ def write_xds_list(xds_list, ref_xds_list, ms_path, output_opts): rechunk=True # Needed to ensure zarr chunks map correctly to disk. ) + # TODO: This needs to be controlled in a sensible way. # write_xds_list = xds_to_table_fragment( # xds_list, # "delta1.ms", diff --git a/quartical/data_handling/predict.py b/quartical/data_handling/predict.py index cc8bda04..e3015e5e 100644 --- a/quartical/data_handling/predict.py +++ b/quartical/data_handling/predict.py @@ -7,7 +7,7 @@ import dask from xarray import DataArray, Dataset from dask.graph_manipulation import clone -from daskms.experimental.multisource import xds_from_table_fragment +from daskms.experimental.fragments import xds_from_table_fragment from loguru import logger import numpy as np From f56b2cf82612550c0d0b7598bd80ca3b483b960f Mon Sep 17 00:00:00 2001 From: landmanbester Date: Thu, 10 Aug 2023 11:26:16 +0200 Subject: [PATCH 24/63] depend on stimela@FIASCO3 and dask-ms@multisource-experimental --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ceee48d2..b6c95c11 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ columnar = "^1.4.1" "ruamel.yaml" = "^0.17.26" dask = {extras = ["diagnostics"], version = "^2023.1.0"} distributed = "^2023.1.0" -dask-ms = {extras = ["s3", "xarray", "zarr"], version = "^0.2.16"} +dask-ms = {git = "https://github.com/ratt-ru/dask-ms.git", branch = "multisource-experimental", extras = ["s3", "xarray", "zarr"]} codex-africanus = {extras = ["dask", "scipy", "astropy", "python-casacore"], version = "^0.3.4"} astro-tigger-lsm = "^1.7.2" loguru = "^0.7.0" @@ -39,7 +39,7 @@ requests = "^2.31.0" pytest = "^7.3.1" omegaconf = "^2.3.0" colorama = "^0.4.6" -stimela = "2.0rc4" +stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO3"} ducc0 = "^0.31.0" sympy = "^1.12" From 037a40d8417453381dd162d9a3ccd1bc7294ca8f Mon Sep 17 00:00:00 2001 From: landmanbester Date: Thu, 10 Aug 2023 12:37:24 +0200 Subject: [PATCH 25/63] depend on stimela@FIASCO3-hlb --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b6c95c11..29e6f275 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,7 +39,7 @@ requests = "^2.31.0" pytest = "^7.3.1" omegaconf = "^2.3.0" colorama = "^0.4.6" -stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO3"} +stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO3-hlb"} ducc0 = "^0.31.0" sympy = "^1.12" From 1bb7448996ec60f273032df48876b4d54d0b2ee4 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Thu, 10 Aug 2023 14:55:40 +0200 Subject: [PATCH 26/63] Add output.fragment_path, which, if set, causes QC to write columns to a dask-ms fragment. --- quartical/config/argument_schema.yaml | 9 ++++++++ quartical/data_handling/ms_handler.py | 30 ++++++++++++++------------- 2 files changed, 25 insertions(+), 14 deletions(-) diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index ac7d768c..7fa416a0 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -193,6 +193,15 @@ output: Name of directory in which QuartiCal logging outputs will be stored. s3 is not currently supported for these outputs. + fragment_path: + dtype: Optional[str] + info: + If set, instead of mutating the input by e.g. writing flags, instead + writes a fragment to this location. A fragment is a zarr backed data + format that is read and dynamically combined with any parent datasets. + This allows QuartiCal to operate in an entirely read-only fashion. + This option is experimental. + log_to_terminal: default: true dtype: bool diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index e1682cfd..4768d62e 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -298,21 +298,23 @@ def write_xds_list(xds_list, ref_xds_list, ms_path, output_opts): with warnings.catch_warnings(): # We anticipate spurious warnings. warnings.simplefilter("ignore") - write_xds_list = xds_to_storage_table( - xds_list, - ms_path, - columns=output_cols, - rechunk=True # Needed to ensure zarr chunks map correctly to disk. - ) - # TODO: This needs to be controlled in a sensible way. - # write_xds_list = xds_to_table_fragment( - # xds_list, - # "delta1.ms", - # ms_path, - # columns=output_cols, - # rechunk=True # Needed to ensure zarr chunks map correctly to disk. - # ) + if output_opts.fragment_path: + write_xds_list = xds_to_table_fragment( + xds_list, + output_opts.fragment_path, + ms_path, + columns=output_cols, + rechunk=True # Ensure zarr chunks map correctly to disk. + ) + + else: + write_xds_list = xds_to_storage_table( + xds_list, + ms_path, + columns=output_cols, + rechunk=True # Ensure zarr chunks map correctly to disk. + ) return write_xds_list From e73450782faedf9234b048713f9356481abb69e8 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 25 Aug 2023 16:03:11 +0200 Subject: [PATCH 27/63] Add plotting functionality (#290) * Fix version drift. * Bump to 0.2.0 * Initial commit of basic plotting functionality. * Change naming convention. * Improve transform argument. * Simplify transform selection. * Add rudimentary time and frequency selection. * Checkpoint ploter changes. Can now handle scans and spws, but is very slow. * More work on plotter - can now plot datasets in parallel. * Some tidying. * Slightly improve plot speed. Dominant cost is still saving the figures. * Commit some minor changes which speed up figure saving. * Lots of tiny fixes. * Tiny cosmetic changes. * Add custom tick formatter so that plots are the same size regardless. * Add matplotlib dependency. * Rework construction of plotting dictionary. Add a few utility functions which will likely be useful in other places in QC. * Rename variable to avoid confusion. * Fix bug affecting recursive grouping. * Avoid copies in grouping code. * Checkpoint work on extending functionality. * Make plotter more powerful. Add colourization option. Begin simplifying interface. * Allow user specification of colourmap. * Add plotsize parameter. --- pyproject.toml | 3 +- quartical/apps/plotter.py | 335 +++++++++++++++++++++++++++++++++ quartical/utils/collections.py | 15 ++ quartical/utils/datasets.py | 36 ++++ 4 files changed, 388 insertions(+), 1 deletion(-) create mode 100644 quartical/apps/plotter.py create mode 100644 quartical/utils/datasets.py diff --git a/pyproject.toml b/pyproject.toml index 68c1f4ee..a67c5df9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,6 +40,7 @@ pytest = "^7.3.1" omegaconf = "^2.3.0" colorama = "^0.4.6" stimela = "2.0rc4" +matplotlib = "^3.5.1" [tool.poetry.scripts] goquartical = 'quartical.executor:execute' @@ -47,7 +48,7 @@ goquartical-config = 'quartical.config.parser:create_user_config' goquartical-backup = 'quartical.apps.backup:backup' goquartical-restore = 'quartical.apps.backup:restore' goquartical-summary = 'quartical.apps.summary:summary' - +goquartical-plot = 'quartical.apps.plotter:plot' [build-system] requires = ["poetry-core"] diff --git a/quartical/apps/plotter.py b/quartical/apps/plotter.py new file mode 100644 index 00000000..35f18977 --- /dev/null +++ b/quartical/apps/plotter.py @@ -0,0 +1,335 @@ +import argparse +import math +import xarray +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import ticker +import matplotlib.cm as cm +from itertools import product, chain +from daskms.experimental.zarr import xds_from_zarr +from daskms.fsspec_store import DaskMSStore +from concurrent.futures import ProcessPoolExecutor, as_completed +from quartical.utils.collections import flatten +from quartical.utils.datasets import recursive_group_by_attr + + +TRANSFORMS = { + "raw": np.array, + "amplitude": np.abs, + "phase": np.angle, + "real": np.real, + "imag": np.imag +} + + +class CustomFormatter(ticker.ScalarFormatter): + + def __init__(self, *args, precision=None, **kwargs): + + super().__init__(*args, *kwargs) + + self.precision = precision + + def _set_format(self): + # set the format string to format all the ticklabels + if len(self.locs) < 2: + # Temporarily augment the locations with the axis end points. + _locs = [*self.locs, *self.axis.get_view_interval()] + else: + _locs = self.locs + locs = (np.asarray(_locs) - self.offset) / 10. ** self.orderOfMagnitude + loc_range = np.ptp(locs) + # Curvilinear coordinates can yield two identical points. + if loc_range == 0: + loc_range = np.max(np.abs(locs)) + # Both points might be zero. + if loc_range == 0: + loc_range = 1 + if len(self.locs) < 2: + # We needed the end points only for the loc_range calculation. + locs = locs[:-2] + loc_range_oom = int(math.floor(math.log10(loc_range))) + # first estimate: + sigfigs = max(0, 3 - loc_range_oom) + # refined estimate: + thresh = 1e-3 * 10 ** loc_range_oom + while sigfigs >= 0: + if np.abs(locs - np.round(locs, decimals=sigfigs)).max() < thresh: + sigfigs -= 1 + else: + break + sigfigs += 1 + sigfigs = self.precision or sigfigs + self.format = f'%{sigfigs + 3}.{sigfigs}f' + if self._usetex or self._useMathText: + self.format = r'$\mathdefault{%s}$' % self.format + + +def cli(): + + parser = argparse.ArgumentParser( + description="Rudimentary plotter for QuartiCal gain solutions." + ) + + parser.add_argument( + "input_path", + type=DaskMSStore, + help="Path to input gains, e.g. path/to/dir/G. Accepts valid s3 urls." + ) + parser.add_argument( + "output_path", + type=DaskMSStore, + help="Path to desired output location." + ) + parser.add_argument( + "--plot-var", + type=str, + default="gains", + help="Name of data variable to plot." + ) + parser.add_argument( + "--flag-var", + type=str, + default="gain_flags", + help="Name of data variable to use as flags." + ) + parser.add_argument( + "--xaxis", + type=str, + default="gain_time", + choices=("gain_time", "gain_freq", "param_time", "param_freq"), + help="Name of coordinate to use for x-axis." + ) + parser.add_argument( + "--transform", + type=str, + default="raw", + choices=list(TRANSFORMS.keys()), + help="Transform to apply to data before plotting." + ) + parser.add_argument( + "--iter-attrs", + type=str, + nargs="+", + default=["FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER"], + help=( + "Attributes (datasets) over which to iterate. Omission will " + "result in concatenation in the omitted axis i.e. omit " + "SCAN_NUMBER to to include all scans in a single plot." + ) + ) + parser.add_argument( + "--iter-axes", + type=str, + nargs="+", + default=["antenna", "direction", "correlation"], + help=( + "Axes over which to iterate when generating plots i.e. produce a " + "plot per unique combination of the specified axes." + ) + ) + parser.add_argument( + "--mean-axis", + type=str, + default=None, + help=( + "If set, will plot a heavier line to indicate the mean of the " + "plotted quantity along this axis." + ) + ) + parser.add_argument( + "--colourize-axis", + type=str, + default=None, + help="Axis to colour by." + ) + parser.add_argument( + "--time-range", + type=float, + nargs=2, + default=[None], + help="Time range to plot." + ) + parser.add_argument( + "--freq-range", + type=float, + nargs=2, + default=[None], + help="Frequency range to plot." + ) + parser.add_argument( + "--nworker", + type=int, + default=1, + help="Number of processes to use while plotting." + ) + parser.add_argument( + "--colourmap", + type=str, + default="plasma", + help=( + "Colourmap to use with --colourize-axis. Supports all matplotlib " + "colourmaps." + ) + ) + parser.add_argument( + "--fig-size", + type=float, + nargs=2, + default=[5, 5], + help="Figure size in inches. Expects two values, width and height." + ) + return parser.parse_args() + + +def to_plot_dict(xdsl, iter_attrs): + + grouped = recursive_group_by_attr(xdsl, iter_attrs) + + return { + k: xarray.combine_by_coords(v, combine_attrs="drop_conflicts") + for k, v in flatten(grouped).items() + } + + +def _plot(group, xds, args): + + xds = xds.compute(scheduler="single-threaded") + + if args.freq_range or args.time_range: + time_ax, freq_ax = xds[args.plot_var].dims[:2] + xds = xds.sel( + { + time_ax: slice(*args.time_range), + freq_ax: slice(*args.freq_range) + } + ) + + dims = xds[args.plot_var].dims # Dimensions of plot quantity. + assert all(map(lambda x: x in dims, args.iter_axes)), ( + f"Some or all of {args.iter_axes} are not present on " + f"{args.plot_var}." + ) + + # Grab the required transform from the dict. + transform = TRANSFORMS[args.transform] + + # NOTE: This mututates the data variables in place. + data = xds[args.plot_var].values + flags = xds[args.flag_var].values + data[np.where(flags)] = np.nan # Set flagged values to nan (not plotted). + xds = xds.drop_vars(args.flag_var) # No more use for flags. + + # Construct list of lists containing axes over which we iterate i.e. + # produce a plot per combination of these values. + iter_axes_itr = [xds[x].values.tolist() for x in args.iter_axes] + + # Figure out axes included in a single plot. + excluded_dims = {*args.iter_axes, args.xaxis} + agg_axes = [d for d in dims if d not in excluded_dims] + agg_axes_itr = [range(xds.sizes[x]) for x in agg_axes] + + # Figure out axes included in a single plot after taking the mean. + excluded_dims = {*args.iter_axes, args.xaxis, args.mean_axis} + mean_agg_axes = [d for d in dims if d not in excluded_dims] + mean_agg_axes_itr = [range(xds.sizes[x]) for x in mean_agg_axes] + + if args.colourize_axis: + n_colour = xds.sizes[args.colourize_axis] + colourmap = cm.get_cmap(args.colourmap) + colours = [colourmap(i / n_colour) for i in range(n_colour)] + else: + n_colour = 2 + colours = ["k", "r"] + + fig, ax = plt.subplots(figsize=args.fig_size) + + for ia in product(*iter_axes_itr): + + sel = {ax: val for ax, val in zip(args.iter_axes, ia)} + + xda = xds.sel(sel)[args.plot_var] + + ax.clear() + + for aa in product(*agg_axes_itr): + + subsel = {ax: val for ax, val in zip(agg_axes, aa)} + pxda = xda.isel(subsel) + + ax.plot( + pxda[args.xaxis].values, + transform(pxda.values), + color=colours[subsel.get(args.colourize_axis, 0)], + linewidth=0.1 + ) + + if args.mean_axis: + + mxda = xda.mean(args.mean_axis) + + for ma in product(*mean_agg_axes_itr): + + subsel = {ax: val for ax, val in zip(mean_agg_axes, ma)} + pxda = mxda.isel(subsel) + + ax.plot( + pxda[args.xaxis].values, + transform(pxda.values), + color=colours[subsel.get(args.colourize_axis, 1)] + ) + + ax.title.set_text("\n".join([f"{k}: {v}" for k, v in sel.items()])) + ax.title.set_fontsize("medium") + ax.set_xlabel(f"{args.xaxis}") + ax.set_ylabel(f"{args.transform}({ia[-1]})") + + formatter = CustomFormatter(precision=2) + formatter.set_scientific(True) + formatter.set_powerlimits((0, 0)) + ax.yaxis.set_major_formatter(formatter) + ax.yaxis.set_major_locator(ticker.LinearLocator(numticks=5)) + + fig_name = "-".join(map(str, chain.from_iterable(sel.items()))) + + root_subdir = f"{xds.NAME}-{args.plot_var}-{args.transform}" + leaf_subdir = "-".join(group) + subdir_path = f"{root_subdir}/{leaf_subdir}" + + args.output_path.makedirs(subdir_path, exist_ok=True) + + fig.savefig( + f"{args.output_path.full_path}/{subdir_path}/{fig_name}.png", + bbox_inches="tight" # SLOW, but slightly unavoidable. + ) + + plt.close() + + +def plot(): + + args = cli() + + non_colourizable_axes = {*args.iter_axes, args.mean_axis, args.xaxis} + if args.colourize_axis and args.colourize_axis in non_colourizable_axes: + raise ValueError(f"Cannot colourize using axis {args.colourize_axis}.") + + # Path to gain location. + gain_path = DaskMSStore("::".join(args.input_path.url.rsplit("/", 1))) + + xdsl = xds_from_zarr(gain_path) + + # Select only the necessary fields for plotting on each dataset. + xdsl = [xds[[args.plot_var, args.flag_var]] for xds in xdsl] + + # Partitioned dictionary of xarray.Datasets. + xdsd = to_plot_dict(xdsl, args.iter_attrs) + + with ProcessPoolExecutor(max_workers=args.nworker) as ppe: + futures = [ppe.submit(_plot, k, xds, args) for k, xds in xdsd.items()] + + for future in as_completed(futures): + try: + future.result() + except Exception as exc: + print(f"Exception raised in process pool: {exc}") diff --git a/quartical/utils/collections.py b/quartical/utils/collections.py index 31da7d34..42199925 100644 --- a/quartical/utils/collections.py +++ b/quartical/utils/collections.py @@ -1,3 +1,4 @@ +from collections.abc import MutableMapping from collections import defaultdict @@ -9,3 +10,17 @@ def freeze_default_dict(ddict): ddict[k] = freeze_default_dict(v) return dict(ddict) + + +def flatten(dictionary, parent_key=()): + """Flatten dictionary. Adapted from https://stackoverflow.com/a/6027615.""" + items = [] + + for key, value in dictionary.items(): + new_key = parent_key + (key,) if parent_key else (key,) + if isinstance(value, MutableMapping): + items.extend(flatten(value, new_key).items()) + else: + items.append((new_key, value)) + + return dict(items) diff --git a/quartical/utils/datasets.py b/quartical/utils/datasets.py new file mode 100644 index 00000000..537bc3a9 --- /dev/null +++ b/quartical/utils/datasets.py @@ -0,0 +1,36 @@ +from collections.abc import Iterable + + +def group_by_attr(xdsl, attr, default="?"): + """Group list of xarray datasets based on value of attribute.""" + + attr_vals = {xds.attrs.get(attr, default) for xds in xdsl} + + return { + f"{attr}_{attr_val}": [ + xds for xds in xdsl if xds.attrs.get(attr, default) == attr_val + ] + for attr_val in attr_vals + } + + +def _recursive_group_by_attr(partition_dict, attrs): + + for k, v in partition_dict.items(): + partition_dict[k] = group_by_attr(v, attrs[0]) + + if len(attrs[1:]): + _recursive_group_by_attr(partition_dict[k], attrs[1:]) + + +def recursive_group_by_attr(xdsl, keys): + + if not isinstance(keys, Iterable): + keys = [keys] + + group_dict = group_by_attr(xdsl, keys[0]) + + if len(keys[1:]): + _recursive_group_by_attr(group_dict, keys[1:]) + + return group_dict From 2e0b9e8c41a99a61d663781c177678437acb30e8 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Wed, 30 Aug 2023 12:16:04 +0200 Subject: [PATCH 28/63] Fix #293 - OOB access caused by `output.subtract_directions` (#294) * Fix version drift. * Bump to 0.2.0 * Fix #293. --- quartical/calibration/calibrate.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/quartical/calibration/calibrate.py b/quartical/calibration/calibrate.py index 83ed022b..dba5f2f2 100644 --- a/quartical/calibration/calibrate.py +++ b/quartical/calibration/calibrate.py @@ -253,6 +253,18 @@ def make_visibility_output( itr = enumerate(zip(data_xds_list, mapping_xds_list)) + if output_opts.subtract_directions: + n_dir = data_xds_list[0].dims['dir'] # Should be the same on all xdss. + requested = set(output_opts.subtract_directions) + valid = set(range(n_dir)) + invalid = requested - valid + if invalid: + raise ValueError( + f"User has specified output.subtract_directions as " + f"{requested} but the following directions are not present " + f"in the model: {invalid}." + ) + for xds_ind, (data_xds, mapping_xds) in itr: data_col = data_xds.DATA.data model_col = data_xds.MODEL_DATA.data From 9a2879466ecf4b9093c8f4d55da9876392e08fe9 Mon Sep 17 00:00:00 2001 From: Landman Bester Date: Wed, 13 Sep 2023 13:14:04 +0200 Subject: [PATCH 29/63] Namedbackups (#296) * Fix version drift. * Bump to 0.2.0 * Add optional label and single field selection to backup app * remove item instead of pop@index * do not .remove() from xds_list * Simplify using some existing functionality. --------- Co-authored-by: JSKenyon Co-authored-by: landmanbester --- quartical/apps/backup.py | 35 +++++++++++++++++++++++----- quartical/data_handling/selection.py | 12 +++++++--- 2 files changed, 38 insertions(+), 9 deletions(-) diff --git a/quartical/apps/backup.py b/quartical/apps/backup.py index fcb44a69..1a15ea8c 100644 --- a/quartical/apps/backup.py +++ b/quartical/apps/backup.py @@ -1,5 +1,6 @@ import argparse from math import prod, ceil +from quartical.data_handling.selection import filter_xds_list from daskms import xds_from_storage_ms, xds_to_storage_table from daskms.experimental.zarr import xds_to_zarr, xds_from_zarr from daskms.fsspec_store import DaskMSStore @@ -10,8 +11,9 @@ def backup(): parser = argparse.ArgumentParser( description='Backup any Measurement Set column to zarr. Backups will ' - 'be labelled automatically using the current datetime, ' - 'the Measurement Set name and the column name.' + 'be labelled using a combination of the passed in label ' + '(defaults to datetime), the Measurement Set name and ' + 'the column name.' ) parser.add_argument( @@ -33,19 +35,34 @@ def backup(): type=str, help='Name of column to be backed up.' ) + parser.add_argument( + '--label', + type=str, + help='An explicit label to include in the backup name. Defaults to ' + 'datetime at which the backup was created. Full name will be ' + 'given by [label]-[msname]-[column].bkp.qc.' + ) parser.add_argument( '--nthread', type=int, default=1, help='Number of threads to use.' ) + parser.add_argument( + '--field-id', + type=int, + help='Field ID to back up.' + ) args = parser.parse_args() ms_name = args.ms_path.full_path.rsplit("/", 1)[1] column_name = args.column_name - timestamp = time.strftime("%Y%m%d-%H%M%S") + if args.label: + label = args.label + else: + label = time.strftime("%Y%m%d-%H%M%S") # This call exists purely to get the relevant shape and dtype info. data_xds_list = xds_from_storage_ms( @@ -55,8 +72,11 @@ def backup(): group_cols=("FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER"), ) + # Use existing functionality. TODO: Improve and expose DDID selection. + xdso = filter_xds_list(data_xds_list, args.field_id) + # Compute appropriate chunks (256MB by default) to keep zarr happy. - chunks = [chunk_by_size(xds[column_name]) for xds in data_xds_list] + chunks = [chunk_by_size(xds[column_name]) for xds in xdso] # Repeat of above call but now with correct chunking information. data_xds_list = xds_from_storage_ms( @@ -67,9 +87,12 @@ def backup(): chunks=chunks ) + # Use existing functionality. TODO: Improve and expose DDID selection. + xdso = filter_xds_list(data_xds_list, args.field_id) + bkp_xds_list = xds_to_zarr( - data_xds_list, - f"{args.zarr_dir.url}::{timestamp}-{ms_name}-{column_name}.bkp.qc", + xdso, + f"{args.zarr_dir.url}::{label}-{ms_name}-{column_name}.bkp.qc", ) dask.compute(bkp_xds_list, num_workers=args.nthread) diff --git a/quartical/data_handling/selection.py b/quartical/data_handling/selection.py index 4f6f9301..bd0f66c3 100644 --- a/quartical/data_handling/selection.py +++ b/quartical/data_handling/selection.py @@ -1,7 +1,13 @@ -def filter_xds_list(xds_list, fields, ddids): +def filter_xds_list(xds_list, fields=[], ddids=[]): - filter_fields = {"FIELD_ID": fields, - "DATA_DESC_ID": ddids} + # If we specify an int, make it a list. Might be worth improving. + fields = [fields] if isinstance(fields, int) else fields + ddids = [ddids] if isinstance(ddids, int) else ddids + + filter_fields = { + "FIELD_ID": fields, + "DATA_DESC_ID": ddids + } for k, v in filter_fields.items(): fil = filter(lambda xds: getattr(xds, k) in v, xds_list) From d8a971969c0774705493ef69141c082b99c4d89b Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Thu, 21 Sep 2023 10:33:12 +0200 Subject: [PATCH 30/63] Depend on dask-ms master. --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index cba975d1..af3e5f49 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ columnar = "^1.4.1" "ruamel.yaml" = "^0.17.26" dask = {extras = ["diagnostics"], version = "^2023.1.0"} distributed = "^2023.1.0" -dask-ms = {git = "https://github.com/ratt-ru/dask-ms.git", branch = "multisource-experimental", extras = ["s3", "xarray", "zarr"]} +dask-ms = {git = "https://github.com/ratt-ru/dask-ms.git", extras = ["s3", "xarray", "zarr"]} codex-africanus = {extras = ["dask", "scipy", "astropy", "python-casacore"], version = "^0.3.4"} astro-tigger-lsm = "^1.7.2" loguru = "^0.7.0" From d7332572718f423c0d1f7891584edd527c19dabe Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Thu, 21 Sep 2023 10:45:29 +0200 Subject: [PATCH 31/63] Add rechunking in fragments code and ROWID to columns in backup app. --- quartical/apps/backup.py | 4 ++-- quartical/data_handling/ms_handler.py | 7 +++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/quartical/apps/backup.py b/quartical/apps/backup.py index 1a15ea8c..594b7d4b 100644 --- a/quartical/apps/backup.py +++ b/quartical/apps/backup.py @@ -67,7 +67,7 @@ def backup(): # This call exists purely to get the relevant shape and dtype info. data_xds_list = xds_from_storage_ms( args.ms_path, - columns=column_name, + columns=(column_name, "ROWID"), index_cols=("TIME",), group_cols=("FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER"), ) @@ -81,7 +81,7 @@ def backup(): # Repeat of above call but now with correct chunking information. data_xds_list = xds_from_storage_ms( args.ms_path, - columns=column_name, + columns=(column_name, "ROWID"), index_cols=("TIME",), group_cols=("FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER"), chunks=chunks diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index 4768d62e..ca652fc8 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -8,6 +8,7 @@ xds_from_table_fragment, xds_to_table_fragment ) +from daskms.experimental.utils import rechunk_by_size from dask.graph_manipulation import clone from loguru import logger from quartical.weights.weights import initialize_weights @@ -300,6 +301,12 @@ def write_xds_list(xds_list, ref_xds_list, ms_path, output_opts): warnings.simplefilter("ignore") if output_opts.fragment_path: + + xds_list = [ + rechunk_by_size(xds, {'corr'}, only_when_needed=True) + for xds in xds_list + ] + write_xds_list = xds_to_table_fragment( xds_list, output_opts.fragment_path, From b6b700e03b428e01793d3d627b900b42af953b91 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 22 Sep 2023 12:59:47 +0200 Subject: [PATCH 32/63] Expose model component help. --- quartical/config/helper.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/quartical/config/helper.py b/quartical/config/helper.py index 00ebecdb..d3df5e7b 100644 --- a/quartical/config/helper.py +++ b/quartical/config/helper.py @@ -14,6 +14,15 @@ "'G.type=complex' or 'B.direction_dependent=True'." ) +ADVANCED_MODEL_MSG = ( + "This section can be safely ignored if 'input_model.advanced_recipe' is " + "not enabled. If enabled, 'input_model.recipe' will behave like " + "'solver.terms' i.e. each component will generate a new config section " + "with an appropriate name e.g. setting 'input_model.recipe=COMP0' " + "will create a configuration section called COMP0 which supports all " + "the options listed here. These can be set using e.g. 'COMP0.type=column'." +) + HELP_MSG = ( "For full help, use 'goquartical help'. For help with a specific " "section, use e.g. 'goquartical help='[section1, section2]''. Help is " @@ -45,6 +54,10 @@ def print_help(HelpConfig, section_names=None): txt = textwrap.fill(GAIN_MSG, width=columns) help_message += f"{Fore.GREEN}{txt:-^{columns}}\n" + if section_name == "model_component": + txt = textwrap.fill(ADVANCED_MODEL_MSG, width=columns) + help_message += f"{Fore.GREEN}{txt:-^{columns}}\n" + for key, value in section.__helpstr__().items(): option = f"{section_name}.{key}" help_message += f"{Fore.MAGENTA}{option:<}\n" @@ -83,7 +96,8 @@ def help(): additional_config = [ oc.from_dotlist( [ - "input_model.advanced_recipe=model_component", + "input_model.recipe=model_component", + "input_model.advanced_recipe=True", "solver.terms=['gain']" ] ) From bb37e9c1c95d3c6affc94bdf2e27840e5912d197 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Tue, 26 Sep 2023 11:09:24 +0200 Subject: [PATCH 33/63] Fix broken identitiy model behaviour caused by addition of degridder. --- quartical/config/preprocess.py | 2 +- quartical/data_handling/model_handler.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/quartical/config/preprocess.py b/quartical/config/preprocess.py index 13dd69be..099686f3 100644 --- a/quartical/config/preprocess.py +++ b/quartical/config/preprocess.py @@ -60,7 +60,7 @@ def transcribe_legacy_recipe(user_recipe): logger.warning( "input_model.recipe was not supplied. Assuming identity model." ) - return IdentityRecipe(Ingredients(set(), set()), dict()) + return IdentityRecipe(Ingredients(set(), set(), set()), dict()) model_columns = set() sky_models = set() diff --git a/quartical/data_handling/model_handler.py b/quartical/data_handling/model_handler.py index 9c237809..87e81373 100644 --- a/quartical/data_handling/model_handler.py +++ b/quartical/data_handling/model_handler.py @@ -212,7 +212,7 @@ def assign_identity_model(data_xds_list): recipe: A modified Recipe object consistent with this case. """ - ingredients = Ingredients({"__IDENT__"}, set()) + ingredients = Ingredients({"__IDENT__"}, set(), set()) instructions = {0: ["__IDENT__"]} recipe = IdentityRecipe(ingredients, instructions) From 9aa54da946e5b885da5df298b22297c5c5bdd000 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Tue, 3 Oct 2023 10:37:43 +0200 Subject: [PATCH 34/63] Selectively disable MAD flagging criteria (#298) * Fix version drift. * Bump to 0.2.0 * Setting MAD threshold to zero will disable flagging on a given statistic. --- quartical/config/argument_schema.yaml | 11 ++++++----- quartical/flagging/flagging_kernels.py | 6 +++--- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index 5e5ad044..be3566fe 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -288,26 +288,27 @@ mad_flags: Multiplicative factor which determines whether or not a chi-squared value is considered to deviate significantly from the median of a given baseline. Values greater than MAD_bl*threshold_bl will be - flagged. + flagged. Set to zero to disable flagging on this statistic. threshold_global: dtype: float - default: 5 + default: 10 info: Multiplicative factor which determines whether or not a chi-squared value is considered to deviate significantly from the median of a given data chunk. Values greater than MAD*threshold_global will be - flagged. + flagged. Set to zero to disable flagging on this statistic. max_deviation: dtype: float - default: 5 + default: 0 info: Multiplicative factor which determines whether or not the MAD estimate on a given baseline is considered to deviate too much from the global MAD estimate. If the MAD estimate over all baselines is X, and the MAD estimate on a specific baseline is X_bl, baselines for which - X_bl > max_deviation*X will be flagged. + X_bl > max_deviation*X will be flagged. Set to zero to disable flagging + on this statistic. solver: terms: diff --git a/quartical/flagging/flagging_kernels.py b/quartical/flagging/flagging_kernels.py index c6a24a8b..2ad152f5 100644 --- a/quartical/flagging/flagging_kernels.py +++ b/quartical/flagging/flagging_kernels.py @@ -110,9 +110,9 @@ def compute_mad_flags( scale_factor = 1.4826 - bl_threshold2 = bl_threshold ** 2 - gbl_threshold2 = gbl_threshold ** 2 - max_deviation2 = max_deviation ** 2 + bl_threshold2 = bl_threshold ** 2 or np.inf + gbl_threshold2 = gbl_threshold ** 2 or np.inf + max_deviation2 = max_deviation ** 2 or np.inf for corr in range(n_corr): From 80428a0427c51388e3fbcb617bf04ee099bc6897 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Thu, 5 Oct 2023 17:56:21 +0200 Subject: [PATCH 35/63] nifty thread safety woes --- quartical/interpolation/interpolants.py | 164 ++++++++++++------------ quartical/interpolation/interpolate.py | 2 + 2 files changed, 87 insertions(+), 79 deletions(-) diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 5151ed77..b2d77b02 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -411,6 +411,8 @@ def smooth_ampphase(gains, fi, to, fo, + ant, + corr, combine_by_time=True, niter=5, nu0=2.0, @@ -419,6 +421,7 @@ def smooth_ampphase(gains, from ducc0.fft import good_size if not combine_by_time: raise NotImplementedError + # weighted sum over time axis ntimei, nchani, nanti, ndiri, ncorri = gains.shape ntimeo = to.size @@ -444,85 +447,88 @@ def smooth_ampphase(gains, print(wgt[mask].max(), wgt[mask].min()) gain[mask] = gain[mask]/wgt[mask] - # set up correlated field model (hardcoding hypers for now) - npix_f = good_size(int(nchani*padding)) - pospace_large = ift.RGSpace([npix_f]) - pospace_small = ift.RGSpace([nchani]) - spfreq_amp = ift.RGSpace(npix_f) - spfreq_phase = ift.RGSpace(npix_f) - - # amplitude model - cfmakera = ift.CorrelatedFieldMaker('amplitude') - cfmakera.add_fluctuations(spfreq_amp, (0.1, 1e-2), None, None, (-3, 1), - 'f') - cfmakera.set_amplitude_total_offset(0., (1e-2, 1e-6)) - cf_amp = cfmakera.finalize() - - normalized_amp = cfmakera.get_normalized_amplitudes() - pspec_freq_amp = normalized_amp[0]**2 - - # phase model - cfmakerp = ift.CorrelatedFieldMaker('phase') - cfmakerp.add_fluctuations(spfreq_phase, (0.1, 1e-2), None, None, (-3, 1), - 'f') - cfmakerp.set_amplitude_total_offset(0., (1e-2, 1e-6)) - cf_phase = cfmakerp.finalize() - - normalized_phase = cfmakerp.get_normalized_amplitudes() - pspec_freq_phase = normalized_phase[0]**2 - - TF = ift.SliceOperator(cf_amp.target, pospace_small.shape) - - signal = TF @ ift.exp(cf_amp.real + 1j*cf_phase.real) - - cf_amp_small = TF @ cf_amp - cf_phase_small = TF @ cf_phase - - tmp = ift.makeField(signal.target, ~mask) - # TODO - this should include averaging - R = ift.MaskOperator(tmp) - dspace = R.target - data = ift.makeField(dspace, gain[mask]) - signal_response = R(signal) - - # Minimization parameters - n_samples = 3 - n_iterations = 5 - ic_sampling = ift.AbsDeltaEnergyController(name='Sampling', deltaE=0.01, iteration_limit=50) - ic_newton = ift.AbsDeltaEnergyController(name='Newton', deltaE=0.01, iteration_limit=15) - minimizer = ift.NewtonCG(ic_newton, enable_logging=True) - - # Set up likelihood energy and information Hamiltonian - W = wgt[mask] - assert (W > 0).all() - covinv = ift.makeField(dspace, W) - Ninv = ift.DiagonalOperator(covinv, sampling_dtype=complex) - inp = (ift.Adder(data, neg=True) @ signal_response) # required to compute the resid for VariableCovGaussianEnergy - BC = ift.ContractionOperator(dspace, None).adjoint # broadcast scalar over everything (None) - weightop = ift.DiagonalOperator(ift.makeField(dspace, W)) - scalarop = ift.GammaOperator(BC.domain, mean=2.0, var=2).ducktape('icov_scalar') - icovop = weightop @ BC @ scalarop.real - inp = inp.ducktape_left('resid') + icovop.ducktape_left('icov') - likelihood_energy = ift.VariableCovarianceGaussianEnergy(R.target, 'resid', 'icov', - data.dtype) @ inp - - samples = ift.optimize_kl(likelihood_energy, - n_iterations, - n_samples, - minimizer, - ic_sampling, - None) # for GeoVI - - # get the mean - signal_mean = samples.average(signal.force).val - - # interpolate to output frequencies - amp = np.interp(fo, fi, np.abs(signal_mean)) - phase = np.interp(fo, fi, np.angle(signal_mean)) - signal_mean = amp*np.exp(1j*phase) - - # broadcast to output times - signal_mean = np.tile(signal_mean[None, :], (to.size, 1)) + with ift.random.Context(int(f'{ant[0]}{corr[0]}')): + + # set up correlated field model (hardcoding hypers for now) + npix_f = good_size(int(nchani*padding)) + pospace_large = ift.RGSpace([npix_f]) + pospace_small = ift.RGSpace([nchani]) + spfreq_amp = ift.RGSpace(npix_f) + spfreq_phase = ift.RGSpace(npix_f) + + # amplitude model + cfmakera = ift.CorrelatedFieldMaker('amplitude') + cfmakera.add_fluctuations(spfreq_amp, (0.1, 1e-2), None, None, (-3, 1), + 'f') + cfmakera.set_amplitude_total_offset(0., (1e-2, 1e-6)) + cf_amp = cfmakera.finalize() + + normalized_amp = cfmakera.get_normalized_amplitudes() + pspec_freq_amp = normalized_amp[0]**2 + + # phase model + cfmakerp = ift.CorrelatedFieldMaker('phase') + cfmakerp.add_fluctuations(spfreq_phase, (0.1, 1e-2), None, None, (-3, 1), + 'f') + cfmakerp.set_amplitude_total_offset(0., (1e-2, 1e-6)) + cf_phase = cfmakerp.finalize() + + normalized_phase = cfmakerp.get_normalized_amplitudes() + pspec_freq_phase = normalized_phase[0]**2 + + TF = ift.SliceOperator(cf_amp.target, pospace_small.shape) + + signal = TF @ ift.exp(cf_amp.real + 1j*cf_phase.real) + + cf_amp_small = TF @ cf_amp + cf_phase_small = TF @ cf_phase + + tmp = ift.makeField(signal.target, ~mask) + # TODO - this should include averaging + R = ift.MaskOperator(tmp) + dspace = R.target + data = ift.makeField(dspace, gain[mask]) + signal_response = R(signal) + + # Minimization parameters + n_samples = 3 + n_iterations = 5 + ic_sampling = ift.AbsDeltaEnergyController(name='Sampling', deltaE=0.01, iteration_limit=50) + ic_newton = ift.AbsDeltaEnergyController(name='Newton', deltaE=0.01, iteration_limit=15) + minimizer = ift.NewtonCG(ic_newton, enable_logging=True) + + # Set up likelihood energy and information Hamiltonian + W = wgt[mask] + assert (W > 0).all() + covinv = ift.makeField(dspace, W) + Ninv = ift.DiagonalOperator(covinv, sampling_dtype=complex) + inp = (ift.Adder(data, neg=True) @ signal_response) # required to compute the resid for VariableCovGaussianEnergy + BC = ift.ContractionOperator(dspace, None).adjoint # broadcast scalar over everything (None) + weightop = ift.DiagonalOperator(ift.makeField(dspace, W)) + scalarop = ift.GammaOperator(BC.domain, mean=2.0, var=2).ducktape('icov_scalar') + icovop = weightop @ BC @ scalarop.real + inp = inp.ducktape_left('resid') + icovop.ducktape_left('icov') + likelihood_energy = ift.VariableCovarianceGaussianEnergy(R.target, 'resid', 'icov', + data.dtype) @ inp + + + samples = ift.optimize_kl(likelihood_energy, + n_iterations, + n_samples, + minimizer, + ic_sampling, + None) # for GeoVI + + # get the mean + signal_mean = samples.average(signal.force).val + + # interpolate to output frequencies + amp = np.interp(fo, fi, np.abs(signal_mean)) + phase = np.interp(fo, fi, np.angle(signal_mean)) + signal_mean = amp*np.exp(1j*phase) + + # broadcast to output times + signal_mean = np.tile(signal_mean[None, :], (to.size, 1)) # broadcast to expected output shape return signal_mean[:, :, None, None, None] diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index d7be4ae9..89501408 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -311,6 +311,8 @@ def bsmooth(merged_xds, target_xds): merged_xds.gain_freq.data, None, out_time, None, out_freq, None, + da.arange(nant, chunks=1), 'p', + da.arange(ncorr, chunks=1), 'c', combine_by_time, None, dtype=merged_xds.gains.dtype, align_arrays=False, From d2808a233c6c0f6e76443d32081f6235ac58099c Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Thu, 19 Oct 2023 13:08:49 +0200 Subject: [PATCH 36/63] Disable mad flagging on off-diagonals by default (#300) * Fix version drift. * Bump to 0.2.0 * Disable flagging based on off-diagonal correlations in the mad flagger by default. This should make the mad flagger less agressive on data with unmodelled polarised emission. --- quartical/config/argument_schema.yaml | 9 +++++++++ quartical/flagging/flagging.py | 10 ++++++++++ quartical/flagging/flagging_kernels.py | 3 ++- 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index be3566fe..52c78af9 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -310,6 +310,15 @@ mad_flags: X_bl > max_deviation*X will be flagged. Set to zero to disable flagging on this statistic. + use_off_diagonals: + dtype: bool + default: False + info: + Controls whether or not the mad flagger will be run on the off-diagonal + elements of the residual. This is disabled by default as the residual + will tend to contain structure in the absence of a polarised model and + adequate leakage calibration. + solver: terms: dtype: List[str] diff --git a/quartical/flagging/flagging.py b/quartical/flagging/flagging.py index 9755c40b..cc4ba1ef 100644 --- a/quartical/flagging/flagging.py +++ b/quartical/flagging/flagging.py @@ -123,6 +123,8 @@ def valid_median(arr): def add_mad_graph(data_xds_list, mad_opts): + diag_corrs = ['RR', 'LL', 'XX', 'YY'] + bl_thresh = mad_opts.threshold_bl gbl_thresh = mad_opts.threshold_global max_deviation = mad_opts.max_deviation @@ -147,6 +149,13 @@ def add_mad_graph(data_xds_list, mad_opts): n_bl_w_autos = (n_ant * (n_ant - 1))/2 + n_ant n_t_chunk, n_f_chunk, _ = residuals.numblocks + if mad_opts.use_off_diagonals: + corr_sel = tuple(np.arange(residuals.shape[-1])) + else: + corr_sel = tuple( + [i for i, c in enumerate(xds.corr.values) if c in diag_corrs] + ) + wres = da.blockwise( compute_whitened_residual, ("rowlike", "chan", "corr"), residuals, ("rowlike", "chan", "corr"), @@ -228,6 +237,7 @@ def add_mad_graph(data_xds_list, mad_opts): gbl_thresh, None, bl_thresh, None, max_deviation, None, + corr_sel, None, n_ant, None, dtype=np.int8, align_arrays=False, diff --git a/quartical/flagging/flagging_kernels.py b/quartical/flagging/flagging_kernels.py index 2ad152f5..21d7eb8a 100644 --- a/quartical/flagging/flagging_kernels.py +++ b/quartical/flagging/flagging_kernels.py @@ -99,6 +99,7 @@ def compute_mad_flags( gbl_threshold, bl_threshold, max_deviation, + corr_sel, n_ant ): @@ -114,7 +115,7 @@ def compute_mad_flags( gbl_threshold2 = gbl_threshold ** 2 or np.inf max_deviation2 = max_deviation ** 2 or np.inf - for corr in range(n_corr): + for corr in corr_sel: gbl_mad_real = gbl_mad_and_med_real[0, 0, corr, 0] * scale_factor gbl_med_real = gbl_mad_and_med_real[0, 0, corr, 1] From 1626f4ab9526d1cebf7600ebf815a227cbe3e186 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Fri, 20 Oct 2023 10:31:04 +0200 Subject: [PATCH 37/63] Fix bug affecting non-standard columns in `input_ms.data_column` (#301) * Fix version drift. * Bump to 0.2.0 * Fix a bug afecting the use of non-standard columns in data column input. --- quartical/data_handling/ms_handler.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index 3b923432..e43d2963 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -89,6 +89,10 @@ def read_xds_list(model_columns, ms_opts): if ms_opts.weight_column not in known_weight_cols: schema[ms_opts.weight_column] = {'dims': ('chan', 'corr')} + known_data_cols = ("DATA", "CORRECTED_DATA") + if ms_opts.data_column not in known_data_cols: + schema[ms_opts.data_column] = {'dims': ('chan', 'corr')} + try: data_xds_list = xds_from_storage_ms( ms_opts.path, From f8de83e52f23ef06968c4d021d9cc3894e482df1 Mon Sep 17 00:00:00 2001 From: Landman Bester Date: Wed, 25 Oct 2023 11:13:53 +0200 Subject: [PATCH 38/63] Don't allow restore app to overwrite metadata (#307) * assign to ms to avoid over-writing metadata in restore app * zip datasets in enumerate * add comment to document failure case * use backup_column_name in restore app * Apply OCD. --------- Co-authored-by: landmanbester Co-authored-by: JSKenyon --- quartical/apps/backup.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/quartical/apps/backup.py b/quartical/apps/backup.py index 1a15ea8c..a5325259 100644 --- a/quartical/apps/backup.py +++ b/quartical/apps/backup.py @@ -135,9 +135,26 @@ def restore(): zarr_root, zarr_name = args.zarr_path.url.rsplit("/", 1) zarr_xds_list = xds_from_zarr(f"{zarr_root}::{zarr_name}") + backup_column_name = list(zarr_xds_list[0].data_vars.keys()).pop() + + # This will fail if the column does not exist but if we allow all columns + # we need to select out the relevant dims for rechunking below + ms_xds_list = xds_from_storage_ms( + args.ms_path, + columns=(args.column_name,), + index_cols=("TIME",), + group_cols=("FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER"), + ) + + for i, (ds, dsr) in enumerate(zip(ms_xds_list, zarr_xds_list)): + dsr = dsr.chunk(ds.chunks) + data_array = getattr(dsr, backup_column_name) + ms_xds_list[i] = ds.assign( + {args.column_name: (data_array.dims, data_array.data)} + ) restored_xds_list = xds_to_storage_table( - zarr_xds_list, + ms_xds_list, args.ms_path, columns=(args.column_name,), rechunk=True From 5a9f7e56339d5442308cbff0a9a1daa502ec7fe9 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 30 Oct 2023 12:25:47 +0200 Subject: [PATCH 39/63] add nifty dependency --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 7d706b4d..982d7b36 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,10 +39,11 @@ requests = "^2.31.0" pytest = "^7.3.1" omegaconf = "^2.3.0" colorama = "^0.4.6" -stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO3-hlb"} +stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO3"} ducc0 = "^0.31.0" sympy = "^1.12" matplotlib = "^3.5.1" +nifty = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8"} [tool.poetry.extras] degrid = ["ducc0", "sympy"] From 0808c00f248fe839b59f6957db2ced0af53a3818 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 30 Oct 2023 12:47:26 +0200 Subject: [PATCH 40/63] nifty -> nifty8 in pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 982d7b36..8e23cae5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,7 +43,7 @@ stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO ducc0 = "^0.31.0" sympy = "^1.12" matplotlib = "^3.5.1" -nifty = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8"} +nifty8 = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8"} [tool.poetry.extras] degrid = ["ducc0", "sympy"] From 8a16e947af0ded1ff15905d0e27b4da32afd969f Mon Sep 17 00:00:00 2001 From: landmanbester Date: Thu, 2 Nov 2023 16:17:44 +0200 Subject: [PATCH 41/63] make concat_by_time optional in 1dsmooth --- quartical/config/argument_schema.yaml | 2 +- quartical/interpolation/interpolants.py | 82 +++++++++++++++---------- quartical/interpolation/interpolate.py | 6 +- 3 files changed, 53 insertions(+), 37 deletions(-) diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index 4ef4c73a..27f8e857 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -438,7 +438,7 @@ dask: scheduler: default: threads - choices: ["threads", "single-threaded", "distributed"] + choices: ["processes", "threads", "single-threaded", "distributed"] info: Which dask scheduler to use. The default, threads, is the most appropriate for non-cluster use. diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index ef5bdd18..c3f57fdf 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -434,43 +434,54 @@ def smooth_ampphase(gains, fo, ant, corr, - combine_by_time=True, + combine_by_time=False, niter=5, nu0=2.0, padding=1.2): + ''' + Use MGVI in nifty to smooth amplitudes and phases. + Here to and fo are the locations where we want to recosntruct the field + whereas ti and fi are the locations where we have data. + ''' import nifty8 as ift from ducc0.fft import good_size - if not combine_by_time: - raise NotImplementedError # weighted sum over time axis ntimei, nchani, nanti, ndiri, ncorri = gains.shape ntimeo = to.size nfreqo = fo.size - assert nanti == 1 + assert ti.size == ntimei + assert fi.size == nchani + assert nanti == 1 # chunked to 1 assert ndiri == 1 - assert ncorri == 1 - gains = gains.squeeze() - jhj = jhj.squeeze().real - gain = np.zeros((nchani), dtype=gains.dtype) - wgt = np.zeros((nchani), dtype=jhj.dtype) - flags = flags.squeeze().astype(bool) + assert ncorri == 1 # chunked to 1 + gains = gains[:, :, 0, 0, 0] + jhj = jhj[:, :, 0, 0, 0].real + flags = flags[:, :, 0, 0].astype(bool) jhj[flags] = 0.0 - for t in range(ntimei): - gain += jhj[t].squeeze() * gains[t].squeeze() - wgt += jhj[t].squeeze() - # normalise by sum of weights - mask = wgt > 0 - if not mask.any(): - return np.ones((ntimeo, nfreqo, nanti, ndiri, ncorri), - dtype=gains.dtype) - - print(wgt[mask].max(), wgt[mask].min()) - gain[mask] = gain[mask]/wgt[mask] - - with ift.random.Context(int(f'{ant[0]}{corr[0]}')): - + if combine_by_time: + # keep time axis + gain = gains[0:1] + wgt = jhj[0:1] + for t in range(1, ntimei): + gain[0:1] += jhj[t:t+1] * gains[t:t+1] + wgt[0:1] += jhj[t:t+1] + # normalise by sum of weights + mask = wgt > 0 + gain[mask] = gain[mask]/wgt[mask] + raise NotImplementedError + else: + gain = gains + wgt = jhj + + ntsmooth = gain.shape[0] + output_gains = np.ones((ntimeo, nfreqo), dtype=gains.dtype) + for t in range(ntsmooth): + mask = wgt[t] > 0 + if not mask.any(): + continue # set up correlated field model (hardcoding hypers for now) + # redo for each t so that each reduction is randomly initialised npix_f = good_size(int(nchani*padding)) pospace_large = ift.RGSpace([npix_f]) pospace_small = ift.RGSpace([nchani]) @@ -505,21 +516,21 @@ def smooth_ampphase(gains, cf_phase_small = TF @ cf_phase tmp = ift.makeField(signal.target, ~mask) - # TODO - this should include averaging + # TODO - this should include down-sampling R = ift.MaskOperator(tmp) dspace = R.target - data = ift.makeField(dspace, gain[mask]) + data = ift.makeField(dspace, gain[t, mask]) signal_response = R(signal) # Minimization parameters n_samples = 3 n_iterations = 5 - ic_sampling = ift.AbsDeltaEnergyController(name='Sampling', deltaE=0.01, iteration_limit=50) - ic_newton = ift.AbsDeltaEnergyController(name='Newton', deltaE=0.01, iteration_limit=15) + ic_sampling = ift.AbsDeltaEnergyController(deltaE=0.01, iteration_limit=50) + ic_newton = ift.AbsDeltaEnergyController(deltaE=0.01, iteration_limit=15) minimizer = ift.NewtonCG(ic_newton, enable_logging=True) # Set up likelihood energy and information Hamiltonian - W = wgt[mask] + W = wgt[t, mask] assert (W > 0).all() covinv = ift.makeField(dspace, W) Ninv = ift.DiagonalOperator(covinv, sampling_dtype=complex) @@ -532,7 +543,6 @@ def smooth_ampphase(gains, likelihood_energy = ift.VariableCovarianceGaussianEnergy(R.target, 'resid', 'icov', data.dtype) @ inp - samples = ift.optimize_kl(likelihood_energy, n_iterations, n_samples, @@ -544,12 +554,18 @@ def smooth_ampphase(gains, signal_mean = samples.average(signal.force).val # interpolate to output frequencies + # TODO - incorporate in response amp = np.interp(fo, fi, np.abs(signal_mean)) phase = np.interp(fo, fi, np.angle(signal_mean)) signal_mean = amp*np.exp(1j*phase) - # broadcast to output times - signal_mean = np.tile(signal_mean[None, :], (to.size, 1)) + if combine_by_time: + # broadcast to output times + output_gains = np.tile(signal_mean[None, :], (to.size, 1)) + else: + output_gains[t] = signal_mean # broadcast to expected output shape - return signal_mean[:, :, None, None, None] + return output_gains[:, :, None, None, None] + + diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index 89501408..4a74e9a4 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -98,7 +98,8 @@ def load_and_interpolate_gains(gain_xds_lod, chain, output_directory): ntime, nchan, nant, ndir, ncorr = merged_xds.gains.shape if ndir > 1: raise ValueError("Smoothing only supported for direction independent gains") - interpolated_xds_list = bsmooth(merged_xds, term_xds_list) + interpolated_xds_list = bsmooth(merged_xds, term_xds_list, + combine_by_time=False) else: # Create a converter object to handle moving between native and # interpolation representations. @@ -285,7 +286,7 @@ def compute_and_reload(directory, gain_xds_list): return gain_xds_list -def bsmooth(merged_xds, target_xds): +def bsmooth(merged_xds, target_xds, combine_by_time=False): out_time = [] out_freq =[] for ds in target_xds: @@ -301,7 +302,6 @@ def bsmooth(merged_xds, target_xds): merged_xds = merged_xds.chunk({"correlation": 1}) from quartical.interpolation.interpolants import smooth_ampphase - combine_by_time = True smoothed_gains = da.blockwise( smooth_ampphase, 'tfpdc', merged_xds.gains.data, 'tfpdc', From e05ceb3e9163ac60f4ed9ef69dd43093059fee82 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Fri, 3 Nov 2023 12:31:19 +0200 Subject: [PATCH 42/63] redirect nifty logs to a file --- quartical/interpolation/interpolants.py | 28 +++++++++++++++++++------ quartical/interpolation/interpolate.py | 6 ++++-- 2 files changed, 26 insertions(+), 8 deletions(-) diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index c3f57fdf..bf3b32e1 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -1,10 +1,15 @@ # -*- coding: utf-8 -*- from loguru import logger # noqa +import logging +import nifty8 as ift +from nifty8.logger import logger +from ducc0.fft import good_size import dask.array as da import numpy as np from scipy.interpolate import RegularGridInterpolator as RGI from numba import njit from quartical.utils.numba import JIT_OPTIONS +from pathlib import Path def linear2d_interpolate_gains(source_xds, target_xds): @@ -434,7 +439,8 @@ def smooth_ampphase(gains, fo, ant, corr, - combine_by_time=False, + combine_by_time, + output_directory, niter=5, nu0=2.0, padding=1.2): @@ -443,8 +449,16 @@ def smooth_ampphase(gains, Here to and fo are the locations where we want to recosntruct the field whereas ti and fi are the locations where we have data. ''' - import nifty8 as ift - from ducc0.fft import good_size + # remove all logger handles and set only a file handler + for handler in logger.handlers: + logger.removeHandler(handler) + logdir = f'{output_directory}/nifty_report_{ant}_{corr}' + path = Path(logdir) + path.mkdir(parents=True, exist_ok=True) + fh = logging.FileHandler(f'{logdir}/output.log', mode='w') + formatter = logging.Formatter('%(asctime)s %(levelname)s %(message)s') + fh.setFormatter(formatter) + logger.addHandler(fh) # weighted sum over time axis ntimei, nchani, nanti, ndiri, ncorri = gains.shape @@ -475,6 +489,8 @@ def smooth_ampphase(gains, wgt = jhj ntsmooth = gain.shape[0] + # TODO - generalise to cover cases for which 1 < ntsmooth < ntimeo + broadcast = ntsmooth == 1 output_gains = np.ones((ntimeo, nfreqo), dtype=gains.dtype) for t in range(ntsmooth): mask = wgt[t] > 0 @@ -559,12 +575,12 @@ def smooth_ampphase(gains, phase = np.interp(fo, fi, np.angle(signal_mean)) signal_mean = amp*np.exp(1j*phase) - if combine_by_time: - # broadcast to output times + if broadcast: output_gains = np.tile(signal_mean[None, :], (to.size, 1)) else: output_gains[t] = signal_mean - + fh.flush() + fh.close() # broadcast to expected output shape return output_gains[:, :, None, None, None] diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index 4a74e9a4..5daac335 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -98,7 +98,7 @@ def load_and_interpolate_gains(gain_xds_lod, chain, output_directory): ntime, nchan, nant, ndir, ncorr = merged_xds.gains.shape if ndir > 1: raise ValueError("Smoothing only supported for direction independent gains") - interpolated_xds_list = bsmooth(merged_xds, term_xds_list, + interpolated_xds_list = bsmooth(merged_xds, term_xds_list, output_directory, combine_by_time=False) else: # Create a converter object to handle moving between native and @@ -286,7 +286,8 @@ def compute_and_reload(directory, gain_xds_list): return gain_xds_list -def bsmooth(merged_xds, target_xds, combine_by_time=False): +def bsmooth(merged_xds, target_xds, output_directory, + combine_by_time=False): out_time = [] out_freq =[] for ds in target_xds: @@ -314,6 +315,7 @@ def bsmooth(merged_xds, target_xds, combine_by_time=False): da.arange(nant, chunks=1), 'p', da.arange(ncorr, chunks=1), 'c', combine_by_time, None, + output_directory, None, dtype=merged_xds.gains.dtype, align_arrays=False, adjust_chunks={'t': ntimeo, 'f': nfreqo}, From 98c51ea49afb4a9075b5b654c5e6ca53f86fe8af Mon Sep 17 00:00:00 2001 From: landmanbester Date: Fri, 3 Nov 2023 12:32:31 +0200 Subject: [PATCH 43/63] propagate dask_opts to interpolation code --- quartical/calibration/calibrate.py | 6 ++++-- quartical/executor.py | 3 ++- quartical/interpolation/interpolate.py | 12 ++++++++---- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/quartical/calibration/calibrate.py b/quartical/calibration/calibrate.py index dba5f2f2..052f68a4 100644 --- a/quartical/calibration/calibrate.py +++ b/quartical/calibration/calibrate.py @@ -114,7 +114,8 @@ def add_calibration_graph( stats_xds_list, solver_opts, chain, - output_opts + output_opts, + dask_opts ): """Given data graph and options, adds the steps necessary for calibration. @@ -159,7 +160,8 @@ def add_calibration_graph( gain_xds_lod = load_and_interpolate_gains( gain_xds_lod, chain, - output_opts.gain_directory + output_opts.gain_directory, + dask_opts ) # Poplulate the gain xarray.Datasets with solutions and convergence info. diff --git a/quartical/executor.py b/quartical/executor.py index b84fb127..fb2b09c9 100644 --- a/quartical/executor.py +++ b/quartical/executor.py @@ -129,7 +129,8 @@ def _execute(exitstack): parangle_xds_list, model_vis_recipe, ms_opts.path, - model_opts + model_opts, + dask_opts ) stats_xds_list = make_stats_xds_list(data_xds_list) diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index 5daac335..c5d893d3 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -16,7 +16,7 @@ from quartical.gains.datasets import write_gain_datasets -def load_and_interpolate_gains(gain_xds_lod, chain, output_directory): +def load_and_interpolate_gains(gain_xds_lod, chain, output_directory, dask_opts): """Load and interpolate gains in accordance with the chain. Given the gain datasets which are to be applied/solved for, determine @@ -133,7 +133,7 @@ def load_and_interpolate_gains(gain_xds_lod, chain, output_directory): # that makes testing painful. interpolated_xds_list = compute_and_reload( - temp_directory, interpolated_xds_list + temp_directory, interpolated_xds_list, dask_opts ) logger.success(f"Successfully loaded/interpolated {term_name}.") @@ -256,7 +256,7 @@ def reindex_and_rechunk(interpolated_xds, reference_xds): return interpolated_xds -def compute_and_reload(directory, gain_xds_list): +def compute_and_reload(directory, gain_xds_list, dask_opts): """Reread gains datasets to be consistent with the reference datasets.""" gain_xds_lod = [{xds.NAME: xds} for xds in gain_xds_list] @@ -264,7 +264,11 @@ def compute_and_reload(directory, gain_xds_list): writes = write_gain_datasets(gain_xds_lod, directory) # NOTE: Need to set compute calls up using dask config mechansim to ensure # correct resource usage is observed. - da.compute(writes) + da.compute(writes, + num_workers=dask_opts.threads, + optimize_graph=True, + scheduler=dask_opts.scheduler, + chunksize=1) # NOTE: This avoids mutating the inputs i.e. avoids side-effects. # TODO: Is this computationally expensive? From 72d8e193095d1b29f46e57398c9c1b8e5f6ecbd6 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Thu, 9 Nov 2023 08:54:57 +0200 Subject: [PATCH 44/63] Fix for summary reporting SOURCE_ID as FIELD_ID (#309) * Fix version drift. * Bump to 0.2.0 * Make summary correctly report FIELD_ID and SOURCE_ID. --- quartical/apps/summary.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/quartical/apps/summary.py b/quartical/apps/summary.py index c6e5d915..36e5e180 100644 --- a/quartical/apps/summary.py +++ b/quartical/apps/summary.py @@ -116,25 +116,32 @@ def field_info(path): field_xds = xds_from_storage_table(path + "::FIELD")[0] - ids = [i for i in field_xds.SOURCE_ID.values] + field_ids = list(range(field_xds.dims['row'])) + source_ids = [i for i in field_xds.SOURCE_ID.values] names = [n for n in field_xds.NAME.values] phase_dirs = [pd for pd in field_xds.PHASE_DIR.values] ref_dirs = [rd for rd in field_xds.REFERENCE_DIR.values] delay_dirs = [dd for dd in field_xds.REFERENCE_DIR.values] msg = "Field summary:\n" - msg += " {:<4} {:<16} {:<16} {:<16} {:<16}\n".format("ID", "NAME", - "PHASE_DIR", - "REF_DIR", - "DELAY_DIR") + msg += " {:<10} {:<10} {:<16} {:<16} {:<16} {:<16}\n".format( + "FIELD_ID", + "SOURCE_ID", + "NAME", + "PHASE_DIR", + "REF_DIR", + "DELAY_DIR" + ) - zipper = zip(ids, names, phase_dirs, ref_dirs, delay_dirs) + zipper = zip( + field_ids, source_ids, names, phase_dirs, ref_dirs, delay_dirs + ) for vals in zipper: - msg += f" {vals[0]:<4} {vals[1]:<16} " \ - f"{'{:.4f} {:.4f}'.format(*vals[2][0]):<16} " \ + msg += f" {vals[0]:<10} {vals[1]:<10} {vals[2]:<16} " \ f"{'{:.4f} {:.4f}'.format(*vals[3][0]):<16} " \ - f"{'{:.4f} {:.4f}'.format(*vals[4][0]):<16}\n" + f"{'{:.4f} {:.4f}'.format(*vals[4][0]):<16} " \ + f"{'{:.4f} {:.4f}'.format(*vals[5][0]):<16}\n" logger.info(msg) From 7b3ee49f7392004ac6c5bd8aa7d754f78cbdd41e Mon Sep 17 00:00:00 2001 From: landmanbester Date: Fri, 10 Nov 2023 12:49:42 +0200 Subject: [PATCH 45/63] move interpolation into executor --- quartical/apps/plotter.py | 12 ++-- quartical/calibration/calibrate.py | 17 +----- quartical/config/argument_schema.yaml | 7 +++ quartical/data_handling/predict.py | 15 +++-- quartical/executor.py | 69 +++++++++++++++++++++-- quartical/interpolation/interpolants.py | 75 +++++++++++++++++-------- quartical/interpolation/interpolate.py | 12 +++- 7 files changed, 151 insertions(+), 56 deletions(-) diff --git a/quartical/apps/plotter.py b/quartical/apps/plotter.py index 35f18977..06149b9b 100644 --- a/quartical/apps/plotter.py +++ b/quartical/apps/plotter.py @@ -216,9 +216,10 @@ def _plot(group, xds, args): # NOTE: This mututates the data variables in place. data = xds[args.plot_var].values - flags = xds[args.flag_var].values - data[np.where(flags)] = np.nan # Set flagged values to nan (not plotted). - xds = xds.drop_vars(args.flag_var) # No more use for flags. + # flags = xds[args.flag_var].values + # data[np.where(flags)] = np.nan # Set flagged values to nan (not plotted). + # xds['flagged_data'] = np.where(flags, data, np.nan) + # xds = xds.drop_vars(args.flag_var) # No more use for flags. # Construct list of lists containing axes over which we iterate i.e. # produce a plot per combination of these values. @@ -244,6 +245,7 @@ def _plot(group, xds, args): fig, ax = plt.subplots(figsize=args.fig_size) + # import ipdb; ipdb.set_trace() for ia in product(*iter_axes_itr): sel = {ax: val for ax, val in zip(args.iter_axes, ia)} @@ -300,6 +302,7 @@ def _plot(group, xds, args): fig.savefig( f"{args.output_path.full_path}/{subdir_path}/{fig_name}.png", + dpi=250, bbox_inches="tight" # SLOW, but slightly unavoidable. ) @@ -320,7 +323,8 @@ def plot(): xdsl = xds_from_zarr(gain_path) # Select only the necessary fields for plotting on each dataset. - xdsl = [xds[[args.plot_var, args.flag_var]] for xds in xdsl] + # xdsl = [xds[[args.plot_var, args.flag_var]] for xds in xdsl] + xdsl = [xds[[args.plot_var]] for xds in xdsl] # Partitioned dictionary of xarray.Datasets. xdsd = to_plot_dict(xdsl, args.iter_attrs) diff --git a/quartical/calibration/calibrate.py b/quartical/calibration/calibrate.py index 052f68a4..aaf9cf6f 100644 --- a/quartical/calibration/calibrate.py +++ b/quartical/calibration/calibrate.py @@ -9,7 +9,6 @@ from quartical.gains.datasets import (make_gain_xds_lod, make_net_xds_lod, populate_net_xds_list) -from quartical.interpolation.interpolate import load_and_interpolate_gains from quartical.gains.baseline import (compute_baseline_corrections, apply_baseline_corrections) from loguru import logger # noqa @@ -114,6 +113,7 @@ def add_calibration_graph( stats_xds_list, solver_opts, chain, + gain_xds_lod, output_opts, dask_opts ): @@ -145,25 +145,10 @@ def add_calibration_graph( "indicate user error." ) - # Create a list of dicts of xarray.Dataset objects which will describe the - # gains per data xarray.Dataset. - gain_xds_lod = make_gain_xds_lod(data_xds_list, chain) - # Create a list of datasets containing mappings. TODO: Is this the best # place to do this? mapping_xds_list = make_mapping_datasets(data_xds_list, chain) - # If there are gains to be loaded from disk, this will load an interpolate - # them to be consistent with this calibration run. TODO: This needs to - # be substantially improved to handle term specific behaviour/utilize - # mappings. - gain_xds_lod = load_and_interpolate_gains( - gain_xds_lod, - chain, - output_opts.gain_directory, - dask_opts - ) - # Poplulate the gain xarray.Datasets with solutions and convergence info. gain_xds_lod, data_xds_list, stats_xds_list = construct_solver( data_xds_list, diff --git a/quartical/config/argument_schema.yaml b/quartical/config/argument_schema.yaml index 27f8e857..14f468e5 100644 --- a/quartical/config/argument_schema.yaml +++ b/quartical/config/argument_schema.yaml @@ -417,6 +417,13 @@ solver: Specify as the integer index of the antenna - antenna names are not currently supported. + bypass: + dtype: bool + default: false + info: + Disable executing the solver graph (eg. if only gain interpolation is + required). + dask: threads: dtype: Optional[int] diff --git a/quartical/data_handling/predict.py b/quartical/data_handling/predict.py index e3015e5e..9c4fb9bb 100644 --- a/quartical/data_handling/predict.py +++ b/quartical/data_handling/predict.py @@ -12,12 +12,15 @@ from loguru import logger import numpy as np import Tigger - -from africanus.util.casa_types import STOKES_ID_MAP -from africanus.util.beams import beam_filenames, beam_grids - -from africanus.experimental.rime.fused import RimeSpecification -from africanus.experimental.rime.fused.dask import rime +from numba.core.errors import NumbaDeprecationWarning +import warnings +with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=NumbaDeprecationWarning) + from africanus.util.casa_types import STOKES_ID_MAP + from africanus.util.beams import beam_filenames, beam_grids + + from africanus.experimental.rime.fused import RimeSpecification + from africanus.experimental.rime.fused.dask import rime from quartical.utils.collections import freeze_default_dict diff --git a/quartical/executor.py b/quartical/executor.py index fb2b09c9..66d089fb 100644 --- a/quartical/executor.py +++ b/quartical/executor.py @@ -19,10 +19,13 @@ from quartical.statistics.logging import log_summary_stats from quartical.flagging.flagging import finalise_flags, add_mad_graph from quartical.scheduling import install_plugin -from quartical.gains.datasets import write_gain_datasets +from quartical.gains.datasets import (write_gain_datasets, + make_gain_xds_lod, + make_net_xds_lod, + populate_net_xds_list) from quartical.gains.baseline import write_baseline_datasets from quartical.utils.dask import compute_context - +from quartical.interpolation.interpolate import load_and_interpolate_gains @logger.catch(onerror=lambda _: sys.exit(1)) def execute(): @@ -129,12 +132,66 @@ def _execute(exitstack): parangle_xds_list, model_vis_recipe, ms_opts.path, - model_opts, - dask_opts + model_opts ) stats_xds_list = make_stats_xds_list(data_xds_list) + # Create a list of dicts of xarray.Dataset objects which will describe the + # gains per data xarray.Dataset. + gain_xds_lod = make_gain_xds_lod(data_xds_list, chain) + + # If there are gains to be loaded from disk, this will load an interpolate + # them to be consistent with this calibration run. TODO: This needs to + # be substantially improved to handle term specific behaviour/utilize + # mappings. + gain_xds_lod = load_and_interpolate_gains( + gain_xds_lod, + chain, + output_opts.gain_directory, + dask_opts + ) + + # if no solving is required then only write out gain datasets. + if solver_opts.bypass: + logger.warning("Bypassing solver!") + + if output_opts.net_gains: + # Construct an effective gain per data_xds. This is always at the full + # time and frequency resolution of the data. Triggers an early compute. + net_xds_lod = make_net_xds_lod( + data_xds_list, + chain, + output_opts + ) + + net_xds_lod = populate_net_xds_list( + net_xds_lod, + gain_xds_lod, + mapping_xds_list, + output_opts + ) + else: + net_xds_lod = [] + + + gain_writes = write_gain_datasets(gain_xds_lod, + output_opts.gain_directory, + net_xds_lod) + + with compute_context(dask_opts, output_opts, time_str): + + dask.compute( + gain_writes, + num_workers=dask_opts.threads, + optimize_graph=True, + scheduler=dask_opts.scheduler + ) + + if dask_opts.scheduler == "distributed": + client.close() # Close this client, hopefully gracefully. + return + # Adds the dask graph describing the calibration of the data. TODO: # This call has excess functionality now. Split out mapping and outputs. cal_outputs = add_calibration_graph( @@ -142,7 +199,9 @@ def _execute(exitstack): stats_xds_list, solver_opts, chain, - output_opts + gain_xds_lod, + output_opts, + dask_opts ) (gain_xds_lod, diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index bf3b32e1..96b67626 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -10,6 +10,7 @@ from numba import njit from quartical.utils.numba import JIT_OPTIONS from pathlib import Path +import warnings def linear2d_interpolate_gains(source_xds, target_xds): @@ -452,7 +453,11 @@ def smooth_ampphase(gains, # remove all logger handles and set only a file handler for handler in logger.handlers: logger.removeHandler(handler) - logdir = f'{output_directory}/nifty_report_{ant}_{corr}' + if isinstance(ant, np.ndarray): + ant=ant[0] + if isinstance(corr, np.ndarray): + corr=corr[0] + logdir = f'{output_directory}/nifty_reports/antenna{ant}_corr{corr}' path = Path(logdir) path.mkdir(parents=True, exist_ok=True) fh = logging.FileHandler(f'{logdir}/output.log', mode='w') @@ -461,11 +466,11 @@ def smooth_ampphase(gains, logger.addHandler(fh) # weighted sum over time axis - ntimei, nchani, nanti, ndiri, ncorri = gains.shape + ntimei, nfreqi, nanti, ndiri, ncorri = gains.shape ntimeo = to.size nfreqo = fo.size assert ti.size == ntimei - assert fi.size == nchani + assert fi.size == nfreqi assert nanti == 1 # chunked to 1 assert ndiri == 1 assert ncorri == 1 # chunked to 1 @@ -474,6 +479,7 @@ def smooth_ampphase(gains, flags = flags[:, :, 0, 0].astype(bool) jhj[flags] = 0.0 if combine_by_time: + # TODO - detrend before combining? # keep time axis gain = gains[0:1] wgt = jhj[0:1] @@ -483,24 +489,21 @@ def smooth_ampphase(gains, # normalise by sum of weights mask = wgt > 0 gain[mask] = gain[mask]/wgt[mask] - raise NotImplementedError else: gain = gains wgt = jhj ntsmooth = gain.shape[0] - # TODO - generalise to cover cases for which 1 < ntsmooth < ntimeo - broadcast = ntsmooth == 1 - output_gains = np.ones((ntimeo, nfreqo), dtype=gains.dtype) + output_gains = np.ones((ntsmooth, nfreqi), dtype=gains.dtype) for t in range(ntsmooth): mask = wgt[t] > 0 if not mask.any(): continue # set up correlated field model (hardcoding hypers for now) # redo for each t so that each reduction is randomly initialised - npix_f = good_size(int(nchani*padding)) + npix_f = good_size(int(nfreqi*padding)) pospace_large = ift.RGSpace([npix_f]) - pospace_small = ift.RGSpace([nchani]) + pospace_small = ift.RGSpace([nfreqi]) spfreq_amp = ift.RGSpace(npix_f) spfreq_phase = ift.RGSpace(npix_f) @@ -532,7 +535,8 @@ def smooth_ampphase(gains, cf_phase_small = TF @ cf_phase tmp = ift.makeField(signal.target, ~mask) - # TODO - this should include down-sampling + # TODO - this should include down-sampling then we can skip the + # frequency interpolation below R = ift.MaskOperator(tmp) dspace = R.target data = ift.makeField(dspace, gain[t, mask]) @@ -559,26 +563,49 @@ def smooth_ampphase(gains, likelihood_energy = ift.VariableCovarianceGaussianEnergy(R.target, 'resid', 'icov', data.dtype) @ inp - samples = ift.optimize_kl(likelihood_energy, - n_iterations, - n_samples, - minimizer, - ic_sampling, - None) # for GeoVI + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=RuntimeWarning) + samples = ift.optimize_kl(likelihood_energy, + n_iterations, + n_samples, + minimizer, + ic_sampling, + None) # for GeoVI # get the mean signal_mean = samples.average(signal.force).val - # interpolate to output frequencies - # TODO - incorporate in response - amp = np.interp(fo, fi, np.abs(signal_mean)) - phase = np.interp(fo, fi, np.angle(signal_mean)) - signal_mean = amp*np.exp(1j*phase) - - if broadcast: + if ntsmooth == 1: + # if there is only a single gain to smooth interpolate + # to output frequencies and broadcast across time + amp = np.interp(fo, fi, np.abs(signal_mean)) + phase = np.interp(fo, fi, np.angle(signal_mean)) + # remove slope and offset + theta = np.polyfit(fo, phase, 1) + phase -= np.polyval(theta, fo) + signal_mean = amp*np.exp(1j*phase) output_gains = np.tile(signal_mean[None, :], (to.size, 1)) else: - output_gains[t] = signal_mean + # otherwise we keep the detrended gains at resolution of the input + # and interpolate afterwards + amp = np.abs(signal_mean) + phase = np.angle(signal_mean) + # remove slope and offset + theta = np.polyfit(fi, phase, 1) + phase -= np.polyval(theta, fi) + output_gains[t] = amp*np.exp(1j*phase) + + # 2D interpolation in this case + # TODO - fill_value=None will extrapolate with linear function but not + # sure this is what we want. Is constant edge exptrapolation possible? + if ntsmooth > 1: + ampo = RGI((ti, fi), np.abs(output_gains), + bounds_error=False, fill_value=None, method='linear') + phaseo = RGI((ti, fi), np.angle(output_gains), + bounds_error=False, fill_value=None, method='linear') + tt, ff = np.meshgrid(to, fo, indexing='ij') + output_gains = ampo((tt, ff)) * np.exp(1j*phaseo((tt, ff))) + fh.flush() fh.close() # broadcast to expected output shape diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index c5d893d3..8f76d1f3 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -264,8 +264,18 @@ def compute_and_reload(directory, gain_xds_list, dask_opts): writes = write_gain_datasets(gain_xds_lod, directory) # NOTE: Need to set compute calls up using dask config mechansim to ensure # correct resource usage is observed. + # import ipdb; ipdb.set_trace() + # from dask import visualize + # name = gain_xds_list[0].NAME + # visualize(writes[0:2], filename=f'/home/bester/projects/ESO137/writes_graph_{name}.pdf', + # optimize_graph=False, engine='cytoscape') + + if dask_opts.scheduler=='distributed': + nworkers = dask_opts.workers + else: + nworkers = dask_opts.threads da.compute(writes, - num_workers=dask_opts.threads, + num_workers=nworkers, optimize_graph=True, scheduler=dask_opts.scheduler, chunksize=1) From 58f9efcaf71b0bb1facbf0cd442583606aff1af6 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Fri, 10 Nov 2023 16:54:33 +0200 Subject: [PATCH 46/63] hack to render params to gains --- quartical/calibration/calibrate.py | 6 +---- quartical/executor.py | 35 +++++++++++++++++++++++++ quartical/interpolation/interpolants.py | 4 ++- 3 files changed, 39 insertions(+), 6 deletions(-) diff --git a/quartical/calibration/calibrate.py b/quartical/calibration/calibrate.py index aaf9cf6f..eeb489b3 100644 --- a/quartical/calibration/calibrate.py +++ b/quartical/calibration/calibrate.py @@ -1,7 +1,6 @@ # -*- coding: utf-8 -*- import numpy as np import dask.array as da -from quartical.calibration.mapping import make_mapping_datasets from quartical.gains.general.generics import (compute_residual, compute_corrected_residual, compute_corrected_weights) @@ -113,6 +112,7 @@ def add_calibration_graph( stats_xds_list, solver_opts, chain, + mapping_xds_list, gain_xds_lod, output_opts, dask_opts @@ -145,10 +145,6 @@ def add_calibration_graph( "indicate user error." ) - # Create a list of datasets containing mappings. TODO: Is this the best - # place to do this? - mapping_xds_list = make_mapping_datasets(data_xds_list, chain) - # Poplulate the gain xarray.Datasets with solutions and convergence info. gain_xds_lod, data_xds_list, stats_xds_list = construct_solver( data_xds_list, diff --git a/quartical/executor.py b/quartical/executor.py index 66d089fb..d5005d5f 100644 --- a/quartical/executor.py +++ b/quartical/executor.py @@ -26,6 +26,7 @@ from quartical.gains.baseline import write_baseline_datasets from quartical.utils.dask import compute_context from quartical.interpolation.interpolate import load_and_interpolate_gains +from quartical.calibration.mapping import make_mapping_datasets @logger.catch(onerror=lambda _: sys.exit(1)) def execute(): @@ -137,6 +138,10 @@ def _execute(exitstack): stats_xds_list = make_stats_xds_list(data_xds_list) + # Create a list of datasets containing mappings. TODO: Is this the best + # place to do this? + mapping_xds_list = make_mapping_datasets(data_xds_list, chain) + # Create a list of dicts of xarray.Dataset objects which will describe the # gains per data xarray.Dataset. gain_xds_lod = make_gain_xds_lod(data_xds_list, chain) @@ -155,6 +160,35 @@ def _execute(exitstack): # if no solving is required then only write out gain datasets. if solver_opts.bypass: logger.warning("Bypassing solver!") + import numpy as np + import dask.array as da + # hack to add gain_flags and convert params to gains + for i, (xdsd, mxds, dxds) in enumerate(zip(gain_xds_lod, + mapping_xds_list, + data_xds_list)): + for key, xds in xdsd.items(): + nt = xds.gain_time.size + nf = xds.gain_freq.size + nant = xds.antenna.size + ndir = xds.direction.size + ncorr = xds.correlation.size + name = xds.NAME + gain_flags = da.zeros((nt, nf, nant, ndir), + chunks=tuple(xds.GAIN_SPEC)[0:-1], + dtype=np.int8) + xds['gain_flags'] = (xds.GAIN_AXES[0:-1], gain_flags) + if xds.TYPE == 'delay': + from quartical.gains.delay.kernel import delay_params_to_gains + params = xds.params.values + gains = np.zeros((nt, nf, nant, ndir, ncorr), np.complex128) + param_freq_map = getattr(mxds, f'{name}_param_freq_map').values + delay_params_to_gains(params, + gains, + dxds.CHAN_FREQ.values, + param_freq_map) + gains = da.from_array(gains, + chunks=tuple(xds.GAIN_SPEC)) + xds['gains'] = (xds.GAIN_AXES, gains) if output_opts.net_gains: # Construct an effective gain per data_xds. This is always at the full @@ -199,6 +233,7 @@ def _execute(exitstack): stats_xds_list, solver_opts, chain, + mapping_xds_list, gain_xds_lod, output_opts, dask_opts diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 96b67626..686e070f 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -1,6 +1,8 @@ # -*- coding: utf-8 -*- from loguru import logger # noqa import logging +import jax.numpy as jnp +logging.getLogger("jax._src.lib.xla_bridge").addFilter(logging.Filter("No GPU/TPU found, falling back to CPU.")) import nifty8 as ift from nifty8.logger import logger from ducc0.fft import good_size @@ -450,7 +452,7 @@ def smooth_ampphase(gains, Here to and fo are the locations where we want to recosntruct the field whereas ti and fi are the locations where we have data. ''' - # remove all logger handles and set only a file handler + # remove all NIFTy logger handles and set only a file handler for handler in logger.handlers: logger.removeHandler(handler) if isinstance(ant, np.ndarray): From c564aa815acf1ffb3cf3c7974a588bd9aec55e18 Mon Sep 17 00:00:00 2001 From: JSKenyon Date: Mon, 13 Nov 2023 08:46:28 +0200 Subject: [PATCH 47/63] Fix receptor summary (#310) * Fix version drift. * Bump to 0.2.0 * Fix incorrect assumption that FEED substable will always have 2 receptors. * Fix similar problem affecting parallactic angle construction. * Update missing column selection for compatibility with upsteam changes. --- quartical/apps/summary.py | 2 +- quartical/data_handling/angles.py | 2 +- quartical/data_handling/ms_handler.py | 28 +++++++++++++++------------ 3 files changed, 18 insertions(+), 14 deletions(-) diff --git a/quartical/apps/summary.py b/quartical/apps/summary.py index 36e5e180..7956a00d 100644 --- a/quartical/apps/summary.py +++ b/quartical/apps/summary.py @@ -99,7 +99,7 @@ def feed_info(path): for i, arrs in enumerate(zipper): for vals in zip(*arrs): msg += f" {i:<4} {vals[0]:<8} {' '.join(vals[1]):<8} " \ - f"{'{:.4f} {:.4f}'.format(*vals[2]):<16}\n" + f"{' '.join([f'{x:.4f}' for x in vals[2]]):<16}\n" logger.info(msg) diff --git a/quartical/data_handling/angles.py b/quartical/data_handling/angles.py index 343ab466..76105d4d 100644 --- a/quartical/data_handling/angles.py +++ b/quartical/data_handling/angles.py @@ -75,7 +75,7 @@ def make_parangle_xds_list(ms_path, data_xds_list): coords={ "utime": np.arange(sum(xds.UTIME_CHUNKS)), "ant": xds.ant, - "receptor": np.arange(2) + "receptor": np.arange(feedtab.dims['receptors']) }, attrs={ "FEED_TYPE": feed_type, diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index e43d2963..fd0e77dd 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -93,19 +93,23 @@ def read_xds_list(model_columns, ms_opts): if ms_opts.data_column not in known_data_cols: schema[ms_opts.data_column] = {'dims': ('chan', 'corr')} - try: - data_xds_list = xds_from_storage_ms( - ms_opts.path, - columns=columns, - index_cols=("TIME",), - group_cols=ms_opts.group_by, - chunks=chunking_per_data_xds, - table_schema=["MS", {**schema}] + data_xds_list = xds_from_storage_ms( + ms_opts.path, + columns=columns, + index_cols=("TIME",), + group_cols=ms_opts.group_by, + chunks=chunking_per_data_xds, + table_schema=["MS", {**schema}] + ) + + missing_columns = set().union( + *[set(columns) - set(xds.data_vars.keys()) for xds in data_xds_list] + ) + + if missing_columns: + raise ValueError( + f"Invalid/missing column specified as input: {missing_columns}." ) - except RuntimeError as e: - raise RuntimeError( - f"Invalid/missing column specified. Underlying error: {e}." - ) from e spw_xds_list = xds_from_storage_table( ms_opts.path + "::SPECTRAL_WINDOW", From 34d2a89f513921a6289e58dbc8a57e4b3bb315ac Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 13 Nov 2023 09:33:36 +0200 Subject: [PATCH 48/63] force branch = 'master' for dask-ms in pyproject.toml --- pyproject.toml | 2 +- quartical/interpolation/interpolants.py | 25 ++++++++++++++----------- quartical/interpolation/interpolate.py | 4 +++- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 8e23cae5..66a19b14 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ columnar = "^1.4.1" "ruamel.yaml" = "^0.17.26" dask = {extras = ["diagnostics"], version = "^2023.1.0"} distributed = "^2023.1.0" -dask-ms = {git = "https://github.com/ratt-ru/dask-ms.git", extras = ["s3", "xarray", "zarr"]} +dask-ms = {git = "https://github.com/ratt-ru/dask-ms.git", extras = ["s3", "xarray", "zarr"], branch="master"} codex-africanus = {extras = ["dask", "scipy", "astropy", "python-casacore"], version = "^0.3.4"} astro-tigger-lsm = "^1.7.2" loguru = "^0.7.0" diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 686e070f..e9747f6e 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -2,7 +2,7 @@ from loguru import logger # noqa import logging import jax.numpy as jnp -logging.getLogger("jax._src.lib.xla_bridge").addFilter(logging.Filter("No GPU/TPU found, falling back to CPU.")) +logging.getLogger("jax._src.xla_bridge").addFilter(logging.Filter("No GPU/TPU found, falling back to CPU.")) import nifty8 as ift from nifty8.logger import logger from ducc0.fft import good_size @@ -577,24 +577,26 @@ def smooth_ampphase(gains, # get the mean signal_mean = samples.average(signal.force).val + amp = np.abs(signal_mean) + phase = np.angle(signal_mean) + # remove slope and offset excluding flagged points + # since the edges could be extrapolated + msk = wgt[t] > 0 + theta = np.polyfit(fi[msk], phase[msk], 1) + print('1') + phase -= np.polyval(theta, fi) + if ntsmooth == 1: # if there is only a single gain to smooth interpolate # to output frequencies and broadcast across time - amp = np.interp(fo, fi, np.abs(signal_mean)) - phase = np.interp(fo, fi, np.angle(signal_mean)) - # remove slope and offset - theta = np.polyfit(fo, phase, 1) - phase -= np.polyval(theta, fo) + amp = np.interp(fo, fi, amp) + phase = np.interp(fo, fi, phase) signal_mean = amp*np.exp(1j*phase) + print('2') output_gains = np.tile(signal_mean[None, :], (to.size, 1)) else: # otherwise we keep the detrended gains at resolution of the input # and interpolate afterwards - amp = np.abs(signal_mean) - phase = np.angle(signal_mean) - # remove slope and offset - theta = np.polyfit(fi, phase, 1) - phase -= np.polyval(theta, fi) output_gains[t] = amp*np.exp(1j*phase) # 2D interpolation in this case @@ -611,6 +613,7 @@ def smooth_ampphase(gains, fh.flush() fh.close() # broadcast to expected output shape + print('3') return output_gains[:, :, None, None, None] diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index 8f76d1f3..599d8df1 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -266,7 +266,9 @@ def compute_and_reload(directory, gain_xds_list, dask_opts): # correct resource usage is observed. # import ipdb; ipdb.set_trace() # from dask import visualize - # name = gain_xds_list[0].NAME + name = gain_xds_list[0].NAME + if name == 'B': + import ipdb; ipdb.set_trace() # visualize(writes[0:2], filename=f'/home/bester/projects/ESO137/writes_graph_{name}.pdf', # optimize_graph=False, engine='cytoscape') From ed922d30a3ca0ce0124a6b7c3e13fc836cb1b753 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 13 Nov 2023 09:44:07 +0200 Subject: [PATCH 49/63] drop python3.8 in ci.yaml --- .github/workflows/ci.yaml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 6590d5c3..24ce4139 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -18,7 +18,7 @@ jobs: strategy: matrix: os: [ubuntu-20.04, ubuntu-22.04] - python-version: ["3.8", "3.9", "3.10"] + python-version: ["3.9", "3.10"] steps: - name: Set up Python ${{ matrix.python-version }} @@ -77,10 +77,10 @@ jobs: # Run on tagged push (generated by tbump). if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags') steps: - - name: Set up Python 3.8 + - name: Set up Python 3.9 uses: actions/setup-python@v4 with: - python-version: 3.8 + python-version: 3.9 - name: Install poetry uses: abatilo/actions-poetry@v2 @@ -120,4 +120,4 @@ jobs: - name: Checkout uses: actions/checkout@v3 - name: Release - uses: softprops/action-gh-release@v0.1.15 \ No newline at end of file + uses: softprops/action-gh-release@v0.1.15 From 14e76d288032fa6d3de95d98947f39456b31cf55 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 13 Nov 2023 09:57:52 +0200 Subject: [PATCH 50/63] fix merge error --- quartical/data_handling/ms_handler.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index c26283c8..2a295c1b 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -106,6 +106,11 @@ def read_xds_list(model_columns, ms_opts): chunks=chunking_per_data_xds, table_schema=["MS", {**schema}] ) + except RuntimeError as e: + raise RuntimeError( + f"Invalid/missing column specified. Underlying error: {e}." + ) from e + spw_xds_list = xds_from_table_fragment( ms_opts.path + "::SPECTRAL_WINDOW", From dfc656312bcee5da2983f66692bf8df5fb035e49 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 13 Nov 2023 11:32:50 +0200 Subject: [PATCH 51/63] python = "^3.9" in pyproject.toml + fix for smoothing gains when full scan is flagged" --- pyproject.toml | 2 +- quartical/interpolation/interpolants.py | 28 ++++++++++--------------- quartical/interpolation/interpolate.py | 11 +++++----- 3 files changed, 18 insertions(+), 23 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 66a19b14..6fd02103 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ include = [ ] [tool.poetry.dependencies] -python = "^3.8" +python = "^3.9" tbump = "^6.10.0" columnar = "^1.4.1" "ruamel.yaml" = "^0.17.26" diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index e9747f6e..8011f5a1 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -496,7 +496,7 @@ def smooth_ampphase(gains, wgt = jhj ntsmooth = gain.shape[0] - output_gains = np.ones((ntsmooth, nfreqi), dtype=gains.dtype) + smooth_gains = np.ones((ntsmooth, nfreqi), dtype=gains.dtype) for t in range(ntsmooth): mask = wgt[t] > 0 if not mask.any(): @@ -583,37 +583,31 @@ def smooth_ampphase(gains, # since the edges could be extrapolated msk = wgt[t] > 0 theta = np.polyfit(fi[msk], phase[msk], 1) - print('1') phase -= np.polyval(theta, fi) - if ntsmooth == 1: - # if there is only a single gain to smooth interpolate - # to output frequencies and broadcast across time - amp = np.interp(fo, fi, amp) - phase = np.interp(fo, fi, phase) - signal_mean = amp*np.exp(1j*phase) - print('2') - output_gains = np.tile(signal_mean[None, :], (to.size, 1)) - else: - # otherwise we keep the detrended gains at resolution of the input - # and interpolate afterwards - output_gains[t] = amp*np.exp(1j*phase) + smooth_gains[t] = amp*np.exp(1j*phase) # 2D interpolation in this case # TODO - fill_value=None will extrapolate with linear function but not # sure this is what we want. Is constant edge exptrapolation possible? if ntsmooth > 1: - ampo = RGI((ti, fi), np.abs(output_gains), + ampo = RGI((ti, fi), np.abs(smooth_gains), bounds_error=False, fill_value=None, method='linear') - phaseo = RGI((ti, fi), np.angle(output_gains), + phaseo = RGI((ti, fi), np.angle(smooth_gains), bounds_error=False, fill_value=None, method='linear') tt, ff = np.meshgrid(to, fo, indexing='ij') output_gains = ampo((tt, ff)) * np.exp(1j*phaseo((tt, ff))) + else: + amp = np.abs(smooth_gains[0]) + amp = np.interp(fo, fi, amp) + phase = np.angle(smooth_gains[0]) + phase = np.interp(fo, fi, phase) + signal_mean = amp*np.exp(1j*phase) + output_gains = np.tile(signal_mean[None, :], (to.size, 1)) fh.flush() fh.close() # broadcast to expected output shape - print('3') return output_gains[:, :, None, None, None] diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index 599d8df1..b0fd28ea 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -266,9 +266,9 @@ def compute_and_reload(directory, gain_xds_list, dask_opts): # correct resource usage is observed. # import ipdb; ipdb.set_trace() # from dask import visualize - name = gain_xds_list[0].NAME - if name == 'B': - import ipdb; ipdb.set_trace() + # name = gain_xds_list[0].NAME + # if name == 'B': + # import ipdb; ipdb.set_trace() # visualize(writes[0:2], filename=f'/home/bester/projects/ESO137/writes_graph_{name}.pdf', # optimize_graph=False, engine='cytoscape') @@ -317,7 +317,7 @@ def bsmooth(merged_xds, target_xds, output_directory, # we also want to chunk by correlation since they can be smoothed separately merged_xds = merged_xds.chunk({"correlation": 1}) - + # import ipdb; ipdb.set_trace() from quartical.interpolation.interpolants import smooth_ampphase smoothed_gains = da.blockwise( smooth_ampphase, 'tfpdc', @@ -344,7 +344,8 @@ def bsmooth(merged_xds, target_xds, output_directory, _, idx0, idx1 = np.intersect1d(t, out_time, assume_unique=True, return_indices=True) - ds['gains'] = (ds.GAIN_AXES, smoothed_gains[idx1]) + target_xds[i] = ds.assign({'gains': (ds.GAIN_AXES, smoothed_gains[idx1])}) + # ds['gains'] = (ds.GAIN_AXES, smoothed_gains[idx1]) target_xds = [ds.chunk({"correlation": ncorr}) for ds in target_xds] From eb257808a7c55b96d3f15d2edfffae8184c94773 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 13 Nov 2023 11:57:13 +0200 Subject: [PATCH 52/63] python = '^3.9, < 3.13' in pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 6fd02103..18df9507 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ include = [ ] [tool.poetry.dependencies] -python = "^3.9" +python = "^3.9, < 3.13" tbump = "^6.10.0" columnar = "^1.4.1" "ruamel.yaml" = "^0.17.26" From 749f271d4bc4f561d0eb2504b542a83bbc0ae8e4 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 13 Nov 2023 12:29:13 +0200 Subject: [PATCH 53/63] add jax nifty dependency --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 18df9507..fa190a78 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,7 +43,7 @@ stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO ducc0 = "^0.31.0" sympy = "^1.12" matplotlib = "^3.5.1" -nifty8 = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8"} +nifty8 = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8", extras = ["re"]} [tool.poetry.extras] degrid = ["ducc0", "sympy"] From 6ac6790758687f3d072d6abc4570d2f7f1dc60bb Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 13 Nov 2023 12:44:40 +0200 Subject: [PATCH 54/63] remove optional jax dependency for nifty as it is not installed correctly --- pyproject.toml | 2 +- quartical/interpolation/interpolants.py | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index fa190a78..18df9507 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,7 +43,7 @@ stimela = {git = "https://github.com/caracal-pipeline/stimela", branch = "FIASCO ducc0 = "^0.31.0" sympy = "^1.12" matplotlib = "^3.5.1" -nifty8 = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8", extras = ["re"]} +nifty8 = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8"} [tool.poetry.extras] degrid = ["ducc0", "sympy"] diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 8011f5a1..97de722f 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -1,8 +1,13 @@ # -*- coding: utf-8 -*- from loguru import logger # noqa import logging -import jax.numpy as jnp -logging.getLogger("jax._src.xla_bridge").addFilter(logging.Filter("No GPU/TPU found, falling back to CPU.")) +# suppress the annoying jax logging if installed +try: + import jax.numpy as jnp + logging.getLogger("jax._src.xla_bridge").addFilter( + logging.Filter("No GPU/TPU found, falling back to CPU.")) +except Exception as e: + pass import nifty8 as ift from nifty8.logger import logger from ducc0.fft import good_size From 5ddd329a6fe972629bea23e04a79522804d3793c Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 14 Nov 2023 10:46:05 +0200 Subject: [PATCH 55/63] pass dask opts in interpolation tests --- quartical/interpolation/interpolants.py | 13 ++++++++++++- testing/tests/interpolation/test_interpolate.py | 6 +++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 97de722f..4de2ae89 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -502,9 +502,19 @@ def smooth_ampphase(gains, ntsmooth = gain.shape[0] smooth_gains = np.ones((ntsmooth, nfreqi), dtype=gains.dtype) - for t in range(ntsmooth): + t = 0 + while t < ntsmooth: mask = wgt[t] > 0 if not mask.any(): + # drop time slot to avoid interpolating to it + # only do this if there is more than one time to smooth + # otherwise there is no array to return + if ntsmooth > 1: + idx = list(np.arange(ntsmooth)) + idx.pop(t) + smooth_gains = smooth_gains[idx] + ntsmooth -= 1 + # do not increment t continue # set up correlated field model (hardcoding hypers for now) # redo for each t so that each reduction is randomly initialised @@ -591,6 +601,7 @@ def smooth_ampphase(gains, phase -= np.polyval(theta, fi) smooth_gains[t] = amp*np.exp(1j*phase) + t += 1 # 2D interpolation in this case # TODO - fill_value=None will extrapolate with linear function but not diff --git a/testing/tests/interpolation/test_interpolate.py b/testing/tests/interpolation/test_interpolate.py index 6948425a..3af38724 100644 --- a/testing/tests/interpolation/test_interpolate.py +++ b/testing/tests/interpolation/test_interpolate.py @@ -112,6 +112,9 @@ def opts(base_opts, interp_mode, interp_method, tmp_path_factory): _opts.G.load_from = str(tmp_path_factory.mktemp("loads.qc")) + "/G" _opts.G.interp_method = interp_method _opts.G.interp_mode = interp_mode + _opts.dask.threads = 8 + _opts.dask.workers = 1 + _opts.dask.scheduler = 'threads' return _opts @@ -181,7 +184,8 @@ def interpolated_xds_lod( return load_and_interpolate_gains( gain_xds_lod, chain, - opts.output.gain_directory + opts.output.gain_directory, + opts.dask ) From 5bb4a093bea45c0efd69591e6e34dab057acf048 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 11 Dec 2023 10:34:23 +0200 Subject: [PATCH 56/63] adjust nifty optimization parameters --- quartical/interpolation/interpolants.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index 4de2ae89..eaaa2101 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -560,10 +560,10 @@ def smooth_ampphase(gains, signal_response = R(signal) # Minimization parameters - n_samples = 3 - n_iterations = 5 - ic_sampling = ift.AbsDeltaEnergyController(deltaE=0.01, iteration_limit=50) - ic_newton = ift.AbsDeltaEnergyController(deltaE=0.01, iteration_limit=15) + n_samples = 10 + n_iterations = 10 + ic_sampling = ift.AbsDeltaEnergyController(deltaE=0.005, iteration_limit=150) + ic_newton = ift.AbsDeltaEnergyController(deltaE=0.005, iteration_limit=25) minimizer = ift.NewtonCG(ic_newton, enable_logging=True) # Set up likelihood energy and information Hamiltonian From ff7ae2929ddf57a513d94ccfcf7c584765a72178 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 19 Dec 2023 11:40:51 +0200 Subject: [PATCH 57/63] make detrending optional during 1dsmooth --- quartical/config/gain_schema.yaml | 6 ++++++ quartical/interpolation/interpolants.py | 10 +++++----- quartical/interpolation/interpolate.py | 5 +++-- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/quartical/config/gain_schema.yaml b/quartical/config/gain_schema.yaml index 070028e2..80c3fb12 100644 --- a/quartical/config/gain_schema.yaml +++ b/quartical/config/gain_schema.yaml @@ -90,3 +90,9 @@ gain: info: Controls whether or not a term will be populated with an initial estimiate where applicable. Currently only supported for delay terms. + + detrend: + dtype: bool + default: true + info: Whether to detrend phases after smoothing. + Only applicable when using interp_method 1dsmooth. diff --git a/quartical/interpolation/interpolants.py b/quartical/interpolation/interpolants.py index eaaa2101..61cd34e3 100644 --- a/quartical/interpolation/interpolants.py +++ b/quartical/interpolation/interpolants.py @@ -449,8 +449,7 @@ def smooth_ampphase(gains, corr, combine_by_time, output_directory, - niter=5, - nu0=2.0, + detrend=True, padding=1.2): ''' Use MGVI in nifty to smooth amplitudes and phases. @@ -596,9 +595,10 @@ def smooth_ampphase(gains, phase = np.angle(signal_mean) # remove slope and offset excluding flagged points # since the edges could be extrapolated - msk = wgt[t] > 0 - theta = np.polyfit(fi[msk], phase[msk], 1) - phase -= np.polyval(theta, fi) + if detrend: + msk = wgt[t] > 0 + theta = np.polyfit(fi[msk], phase[msk], 1) + phase -= np.polyval(theta, fi) smooth_gains[t] = amp*np.exp(1j*phase) t += 1 diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index b0fd28ea..30123696 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -99,7 +99,7 @@ def load_and_interpolate_gains(gain_xds_lod, chain, output_directory, dask_opts) if ndir > 1: raise ValueError("Smoothing only supported for direction independent gains") interpolated_xds_list = bsmooth(merged_xds, term_xds_list, output_directory, - combine_by_time=False) + combine_by_time=False, detrend=term.detrend) else: # Create a converter object to handle moving between native and # interpolation representations. @@ -303,7 +303,7 @@ def compute_and_reload(directory, gain_xds_list, dask_opts): def bsmooth(merged_xds, target_xds, output_directory, - combine_by_time=False): + combine_by_time=False, detrend=True): out_time = [] out_freq =[] for ds in target_xds: @@ -332,6 +332,7 @@ def bsmooth(merged_xds, target_xds, output_directory, da.arange(ncorr, chunks=1), 'c', combine_by_time, None, output_directory, None, + detrend, None, dtype=merged_xds.gains.dtype, align_arrays=False, adjust_chunks={'t': ntimeo, 'f': nfreqo}, From 13d26ccd12eb9e348b7c92b1c155bf4a6c2ee575 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 19 Dec 2023 11:45:12 +0200 Subject: [PATCH 58/63] allow smoothing complex type gains --- quartical/interpolation/interpolate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/quartical/interpolation/interpolate.py b/quartical/interpolation/interpolate.py index 30123696..9a70c7dd 100644 --- a/quartical/interpolation/interpolate.py +++ b/quartical/interpolation/interpolate.py @@ -93,8 +93,8 @@ def load_and_interpolate_gains(gain_xds_lod, chain, output_directory, dask_opts) merged_xds = merged_xds.chunk({**merged_xds.dims, "antenna": 1}) if term.interp_method == "1dsmooth": - if merged_xds.TYPE != "diag_complex": - raise ValueError("Smoothing only supported for diag_complex type") + if merged_xds.TYPE not in ["diag_complex", "complex"]: + raise ValueError(f"Smoothing not suported for term type {merged_xds.TYPE}") ntime, nchan, nant, ndir, ncorr = merged_xds.gains.shape if ndir > 1: raise ValueError("Smoothing only supported for direction independent gains") From dfb390e7f36242ca84cfd662215eb0a4d4cc2ffa Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 30 Jan 2024 11:44:15 +0200 Subject: [PATCH 59/63] python >= 3.10 for nifty_8 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index f248bf22..5a350c6b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ include = [ ] [tool.poetry.dependencies] -python = "^3.9" +python = "^3.10" astro-tigger-lsm = ">=1.7.2, <=1.7.3" codex-africanus = {extras = ["dask", "scipy", "astropy", "python-casacore"], version = ">=0.3.4, <=0.3.4"} colorama = ">=0.4.6, <=0.4.6" From 7379267e387aab4e8333f62416ded7ec2e451821 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 30 Jan 2024 11:49:47 +0200 Subject: [PATCH 60/63] remove 3.9 from testing matrix --- .github/workflows/ci.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 24ce4139..c8078056 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -18,7 +18,7 @@ jobs: strategy: matrix: os: [ubuntu-20.04, ubuntu-22.04] - python-version: ["3.9", "3.10"] + python-version: ["3.10"] steps: - name: Set up Python ${{ matrix.python-version }} @@ -77,10 +77,10 @@ jobs: # Run on tagged push (generated by tbump). if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags') steps: - - name: Set up Python 3.9 + - name: Set up Python 3.10 uses: actions/setup-python@v4 with: - python-version: 3.9 + python-version: 3.10 - name: Install poetry uses: abatilo/actions-poetry@v2 From 23491ab0d7eb32b7aace1d3f1d6b8f6cab2ad056 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 30 Jan 2024 14:28:00 +0200 Subject: [PATCH 61/63] regenerate lock file --- poetry.lock | 298 +++++++++++++++++++++++++++++-------------------- pyproject.toml | 2 + 2 files changed, 179 insertions(+), 121 deletions(-) diff --git a/poetry.lock b/poetry.lock index f9cb2208..3cfcab9c 100644 --- a/poetry.lock +++ b/poetry.lock @@ -23,87 +23,87 @@ boto3 = ["boto3 (>=1.33.2,<1.34.28)"] [[package]] name = "aiohttp" -version = "3.9.2" +version = "3.9.3" description = "Async http client/server framework (asyncio)" optional = false python-versions = ">=3.8" files = [ - {file = "aiohttp-3.9.2-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:772fbe371788e61c58d6d3d904268e48a594ba866804d08c995ad71b144f94cb"}, - {file = "aiohttp-3.9.2-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:edd4f1af2253f227ae311ab3d403d0c506c9b4410c7fc8d9573dec6d9740369f"}, - {file = "aiohttp-3.9.2-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:cfee9287778399fdef6f8a11c9e425e1cb13cc9920fd3a3df8f122500978292b"}, - {file = "aiohttp-3.9.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:3cc158466f6a980a6095ee55174d1de5730ad7dec251be655d9a6a9dd7ea1ff9"}, - {file = "aiohttp-3.9.2-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:54ec82f45d57c9a65a1ead3953b51c704f9587440e6682f689da97f3e8defa35"}, - {file = "aiohttp-3.9.2-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:abeb813a18eb387f0d835ef51f88568540ad0325807a77a6e501fed4610f864e"}, - {file = "aiohttp-3.9.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:cc91d07280d7d169f3a0f9179d8babd0ee05c79d4d891447629ff0d7d8089ec2"}, - {file = "aiohttp-3.9.2-cp310-cp310-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:b65e861f4bebfb660f7f0f40fa3eb9f2ab9af10647d05dac824390e7af8f75b7"}, - {file = "aiohttp-3.9.2-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:04fd8ffd2be73d42bcf55fd78cde7958eeee6d4d8f73c3846b7cba491ecdb570"}, - {file = "aiohttp-3.9.2-cp310-cp310-musllinux_1_1_i686.whl", hash = "sha256:3d8d962b439a859b3ded9a1e111a4615357b01620a546bc601f25b0211f2da81"}, - {file = "aiohttp-3.9.2-cp310-cp310-musllinux_1_1_ppc64le.whl", hash = "sha256:8ceb658afd12b27552597cf9a65d9807d58aef45adbb58616cdd5ad4c258c39e"}, - {file = "aiohttp-3.9.2-cp310-cp310-musllinux_1_1_s390x.whl", hash = "sha256:0e4ee4df741670560b1bc393672035418bf9063718fee05e1796bf867e995fad"}, - {file = "aiohttp-3.9.2-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:2dec87a556f300d3211decf018bfd263424f0690fcca00de94a837949fbcea02"}, - {file = "aiohttp-3.9.2-cp310-cp310-win32.whl", hash = "sha256:3e1a800f988ce7c4917f34096f81585a73dbf65b5c39618b37926b1238cf9bc4"}, - {file = "aiohttp-3.9.2-cp310-cp310-win_amd64.whl", hash = "sha256:ea510718a41b95c236c992b89fdfc3d04cc7ca60281f93aaada497c2b4e05c46"}, - {file = "aiohttp-3.9.2-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:6aaa6f99256dd1b5756a50891a20f0d252bd7bdb0854c5d440edab4495c9f973"}, - {file = "aiohttp-3.9.2-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:a27d8c70ad87bcfce2e97488652075a9bdd5b70093f50b10ae051dfe5e6baf37"}, - {file = "aiohttp-3.9.2-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:54287bcb74d21715ac8382e9de146d9442b5f133d9babb7e5d9e453faadd005e"}, - {file = "aiohttp-3.9.2-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5bb3d05569aa83011fcb346b5266e00b04180105fcacc63743fc2e4a1862a891"}, - {file = "aiohttp-3.9.2-cp311-cp311-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:c8534e7d69bb8e8d134fe2be9890d1b863518582f30c9874ed7ed12e48abe3c4"}, - {file = "aiohttp-3.9.2-cp311-cp311-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:4bd9d5b989d57b41e4ff56ab250c5ddf259f32db17159cce630fd543376bd96b"}, - {file = "aiohttp-3.9.2-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:fa6904088e6642609981f919ba775838ebf7df7fe64998b1a954fb411ffb4663"}, - {file = "aiohttp-3.9.2-cp311-cp311-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:bda42eb410be91b349fb4ee3a23a30ee301c391e503996a638d05659d76ea4c2"}, - {file = "aiohttp-3.9.2-cp311-cp311-musllinux_1_1_aarch64.whl", hash = "sha256:193cc1ccd69d819562cc7f345c815a6fc51d223b2ef22f23c1a0f67a88de9a72"}, - {file = "aiohttp-3.9.2-cp311-cp311-musllinux_1_1_i686.whl", hash = "sha256:b9f1cb839b621f84a5b006848e336cf1496688059d2408e617af33e3470ba204"}, - {file = "aiohttp-3.9.2-cp311-cp311-musllinux_1_1_ppc64le.whl", hash = "sha256:d22a0931848b8c7a023c695fa2057c6aaac19085f257d48baa24455e67df97ec"}, - {file = "aiohttp-3.9.2-cp311-cp311-musllinux_1_1_s390x.whl", hash = "sha256:4112d8ba61fbd0abd5d43a9cb312214565b446d926e282a6d7da3f5a5aa71d36"}, - {file = "aiohttp-3.9.2-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:c4ad4241b52bb2eb7a4d2bde060d31c2b255b8c6597dd8deac2f039168d14fd7"}, - {file = "aiohttp-3.9.2-cp311-cp311-win32.whl", hash = "sha256:ee2661a3f5b529f4fc8a8ffee9f736ae054adfb353a0d2f78218be90617194b3"}, - {file = "aiohttp-3.9.2-cp311-cp311-win_amd64.whl", hash = "sha256:4deae2c165a5db1ed97df2868ef31ca3cc999988812e82386d22937d9d6fed52"}, - {file = "aiohttp-3.9.2-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:6f4cdba12539215aaecf3c310ce9d067b0081a0795dd8a8805fdb67a65c0572a"}, - {file = "aiohttp-3.9.2-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:84e843b33d5460a5c501c05539809ff3aee07436296ff9fbc4d327e32aa3a326"}, - {file = "aiohttp-3.9.2-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:8008d0f451d66140a5aa1c17e3eedc9d56e14207568cd42072c9d6b92bf19b52"}, - {file = "aiohttp-3.9.2-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:61c47ab8ef629793c086378b1df93d18438612d3ed60dca76c3422f4fbafa792"}, - {file = "aiohttp-3.9.2-cp312-cp312-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:bc71f748e12284312f140eaa6599a520389273174b42c345d13c7e07792f4f57"}, - {file = "aiohttp-3.9.2-cp312-cp312-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:a1c3a4d0ab2f75f22ec80bca62385db2e8810ee12efa8c9e92efea45c1849133"}, - {file = "aiohttp-3.9.2-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:9a87aa0b13bbee025faa59fa58861303c2b064b9855d4c0e45ec70182bbeba1b"}, - {file = "aiohttp-3.9.2-cp312-cp312-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:e2cc0d04688b9f4a7854c56c18aa7af9e5b0a87a28f934e2e596ba7e14783192"}, - {file = "aiohttp-3.9.2-cp312-cp312-musllinux_1_1_aarch64.whl", hash = "sha256:1956e3ac376b1711c1533266dec4efd485f821d84c13ce1217d53e42c9e65f08"}, - {file = "aiohttp-3.9.2-cp312-cp312-musllinux_1_1_i686.whl", hash = "sha256:114da29f39eccd71b93a0fcacff178749a5c3559009b4a4498c2c173a6d74dff"}, - {file = "aiohttp-3.9.2-cp312-cp312-musllinux_1_1_ppc64le.whl", hash = "sha256:3f17999ae3927d8a9a823a1283b201344a0627272f92d4f3e3a4efe276972fe8"}, - {file = "aiohttp-3.9.2-cp312-cp312-musllinux_1_1_s390x.whl", hash = "sha256:f31df6a32217a34ae2f813b152a6f348154f948c83213b690e59d9e84020925c"}, - {file = "aiohttp-3.9.2-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:7a75307ffe31329928a8d47eae0692192327c599113d41b278d4c12b54e1bd11"}, - {file = "aiohttp-3.9.2-cp312-cp312-win32.whl", hash = "sha256:972b63d589ff8f305463593050a31b5ce91638918da38139b9d8deaba9e0fed7"}, - {file = "aiohttp-3.9.2-cp312-cp312-win_amd64.whl", hash = "sha256:200dc0246f0cb5405c80d18ac905c8350179c063ea1587580e3335bfc243ba6a"}, - {file = "aiohttp-3.9.2-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:158564d0d1020e0d3fe919a81d97aadad35171e13e7b425b244ad4337fc6793a"}, - {file = "aiohttp-3.9.2-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:da1346cd0ccb395f0ed16b113ebb626fa43b7b07fd7344fce33e7a4f04a8897a"}, - {file = "aiohttp-3.9.2-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:eaa9256de26ea0334ffa25f1913ae15a51e35c529a1ed9af8e6286dd44312554"}, - {file = "aiohttp-3.9.2-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:1543e7fb00214fb4ccead42e6a7d86f3bb7c34751ec7c605cca7388e525fd0b4"}, - {file = "aiohttp-3.9.2-cp38-cp38-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:186e94570433a004e05f31f632726ae0f2c9dee4762a9ce915769ce9c0a23d89"}, - {file = "aiohttp-3.9.2-cp38-cp38-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:d52d20832ac1560f4510d68e7ba8befbc801a2b77df12bd0cd2bcf3b049e52a4"}, - {file = "aiohttp-3.9.2-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1c45e4e815ac6af3b72ca2bde9b608d2571737bb1e2d42299fc1ffdf60f6f9a1"}, - {file = "aiohttp-3.9.2-cp38-cp38-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:aa906b9bdfd4a7972dd0628dbbd6413d2062df5b431194486a78f0d2ae87bd55"}, - {file = "aiohttp-3.9.2-cp38-cp38-musllinux_1_1_aarch64.whl", hash = "sha256:68bbee9e17d66f17bb0010aa15a22c6eb28583edcc8b3212e2b8e3f77f3ebe2a"}, - {file = "aiohttp-3.9.2-cp38-cp38-musllinux_1_1_i686.whl", hash = "sha256:4c189b64bd6d9a403a1a3f86a3ab3acbc3dc41a68f73a268a4f683f89a4dec1f"}, - {file = "aiohttp-3.9.2-cp38-cp38-musllinux_1_1_ppc64le.whl", hash = "sha256:8a7876f794523123bca6d44bfecd89c9fec9ec897a25f3dd202ee7fc5c6525b7"}, - {file = "aiohttp-3.9.2-cp38-cp38-musllinux_1_1_s390x.whl", hash = "sha256:d23fba734e3dd7b1d679b9473129cd52e4ec0e65a4512b488981a56420e708db"}, - {file = "aiohttp-3.9.2-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:b141753be581fab842a25cb319f79536d19c2a51995d7d8b29ee290169868eab"}, - {file = "aiohttp-3.9.2-cp38-cp38-win32.whl", hash = "sha256:103daf41ff3b53ba6fa09ad410793e2e76c9d0269151812e5aba4b9dd674a7e8"}, - {file = "aiohttp-3.9.2-cp38-cp38-win_amd64.whl", hash = "sha256:328918a6c2835861ff7afa8c6d2c70c35fdaf996205d5932351bdd952f33fa2f"}, - {file = "aiohttp-3.9.2-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:5264d7327c9464786f74e4ec9342afbbb6ee70dfbb2ec9e3dfce7a54c8043aa3"}, - {file = "aiohttp-3.9.2-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:07205ae0015e05c78b3288c1517afa000823a678a41594b3fdc870878d645305"}, - {file = "aiohttp-3.9.2-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:ae0a1e638cffc3ec4d4784b8b4fd1cf28968febc4bd2718ffa25b99b96a741bd"}, - {file = "aiohttp-3.9.2-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:d43302a30ba1166325974858e6ef31727a23bdd12db40e725bec0f759abce505"}, - {file = "aiohttp-3.9.2-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:16a967685907003765855999af11a79b24e70b34dc710f77a38d21cd9fc4f5fe"}, - {file = "aiohttp-3.9.2-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:6fa3ee92cd441d5c2d07ca88d7a9cef50f7ec975f0117cd0c62018022a184308"}, - {file = "aiohttp-3.9.2-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:0b500c5ad9c07639d48615a770f49618130e61be36608fc9bc2d9bae31732b8f"}, - {file = "aiohttp-3.9.2-cp39-cp39-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:c07327b368745b1ce2393ae9e1aafed7073d9199e1dcba14e035cc646c7941bf"}, - {file = "aiohttp-3.9.2-cp39-cp39-musllinux_1_1_aarch64.whl", hash = "sha256:cc7d6502c23a0ec109687bf31909b3fb7b196faf198f8cff68c81b49eb316ea9"}, - {file = "aiohttp-3.9.2-cp39-cp39-musllinux_1_1_i686.whl", hash = "sha256:07be2be7071723c3509ab5c08108d3a74f2181d4964e869f2504aaab68f8d3e8"}, - {file = "aiohttp-3.9.2-cp39-cp39-musllinux_1_1_ppc64le.whl", hash = "sha256:122468f6fee5fcbe67cb07014a08c195b3d4c41ff71e7b5160a7bcc41d585a5f"}, - {file = "aiohttp-3.9.2-cp39-cp39-musllinux_1_1_s390x.whl", hash = "sha256:00a9abcea793c81e7f8778ca195a1714a64f6d7436c4c0bb168ad2a212627000"}, - {file = "aiohttp-3.9.2-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:7a9825fdd64ecac5c670234d80bb52bdcaa4139d1f839165f548208b3779c6c6"}, - {file = "aiohttp-3.9.2-cp39-cp39-win32.whl", hash = "sha256:5422cd9a4a00f24c7244e1b15aa9b87935c85fb6a00c8ac9b2527b38627a9211"}, - {file = "aiohttp-3.9.2-cp39-cp39-win_amd64.whl", hash = "sha256:7d579dcd5d82a86a46f725458418458fa43686f6a7b252f2966d359033ffc8ab"}, - {file = "aiohttp-3.9.2.tar.gz", hash = "sha256:b0ad0a5e86ce73f5368a164c10ada10504bf91869c05ab75d982c6048217fbf7"}, + {file = "aiohttp-3.9.3-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:939677b61f9d72a4fa2a042a5eee2a99a24001a67c13da113b2e30396567db54"}, + {file = "aiohttp-3.9.3-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:1f5cd333fcf7590a18334c90f8c9147c837a6ec8a178e88d90a9b96ea03194cc"}, + {file = "aiohttp-3.9.3-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:82e6aa28dd46374f72093eda8bcd142f7771ee1eb9d1e223ff0fa7177a96b4a5"}, + {file = "aiohttp-3.9.3-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f56455b0c2c7cc3b0c584815264461d07b177f903a04481dfc33e08a89f0c26b"}, + {file = "aiohttp-3.9.3-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:bca77a198bb6e69795ef2f09a5f4c12758487f83f33d63acde5f0d4919815768"}, + {file = "aiohttp-3.9.3-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:e083c285857b78ee21a96ba1eb1b5339733c3563f72980728ca2b08b53826ca5"}, + {file = "aiohttp-3.9.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ab40e6251c3873d86ea9b30a1ac6d7478c09277b32e14745d0d3c6e76e3c7e29"}, + {file = "aiohttp-3.9.3-cp310-cp310-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:df822ee7feaaeffb99c1a9e5e608800bd8eda6e5f18f5cfb0dc7eeb2eaa6bbec"}, + {file = "aiohttp-3.9.3-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:acef0899fea7492145d2bbaaaec7b345c87753168589cc7faf0afec9afe9b747"}, + {file = "aiohttp-3.9.3-cp310-cp310-musllinux_1_1_i686.whl", hash = "sha256:cd73265a9e5ea618014802ab01babf1940cecb90c9762d8b9e7d2cc1e1969ec6"}, + {file = "aiohttp-3.9.3-cp310-cp310-musllinux_1_1_ppc64le.whl", hash = "sha256:a78ed8a53a1221393d9637c01870248a6f4ea5b214a59a92a36f18151739452c"}, + {file = "aiohttp-3.9.3-cp310-cp310-musllinux_1_1_s390x.whl", hash = "sha256:6b0e029353361f1746bac2e4cc19b32f972ec03f0f943b390c4ab3371840aabf"}, + {file = "aiohttp-3.9.3-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:7cf5c9458e1e90e3c390c2639f1017a0379a99a94fdfad3a1fd966a2874bba52"}, + {file = "aiohttp-3.9.3-cp310-cp310-win32.whl", hash = "sha256:3e59c23c52765951b69ec45ddbbc9403a8761ee6f57253250c6e1536cacc758b"}, + {file = "aiohttp-3.9.3-cp310-cp310-win_amd64.whl", hash = "sha256:055ce4f74b82551678291473f66dc9fb9048a50d8324278751926ff0ae7715e5"}, + {file = "aiohttp-3.9.3-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:6b88f9386ff1ad91ace19d2a1c0225896e28815ee09fc6a8932fded8cda97c3d"}, + {file = "aiohttp-3.9.3-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:c46956ed82961e31557b6857a5ca153c67e5476972e5f7190015018760938da2"}, + {file = "aiohttp-3.9.3-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:07b837ef0d2f252f96009e9b8435ec1fef68ef8b1461933253d318748ec1acdc"}, + {file = "aiohttp-3.9.3-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:dad46e6f620574b3b4801c68255492e0159d1712271cc99d8bdf35f2043ec266"}, + {file = "aiohttp-3.9.3-cp311-cp311-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:5ed3e046ea7b14938112ccd53d91c1539af3e6679b222f9469981e3dac7ba1ce"}, + {file = "aiohttp-3.9.3-cp311-cp311-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:039df344b45ae0b34ac885ab5b53940b174530d4dd8a14ed8b0e2155b9dddccb"}, + {file = "aiohttp-3.9.3-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:7943c414d3a8d9235f5f15c22ace69787c140c80b718dcd57caaade95f7cd93b"}, + {file = "aiohttp-3.9.3-cp311-cp311-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:84871a243359bb42c12728f04d181a389718710129b36b6aad0fc4655a7647d4"}, + {file = "aiohttp-3.9.3-cp311-cp311-musllinux_1_1_aarch64.whl", hash = "sha256:5eafe2c065df5401ba06821b9a054d9cb2848867f3c59801b5d07a0be3a380ae"}, + {file = "aiohttp-3.9.3-cp311-cp311-musllinux_1_1_i686.whl", hash = "sha256:9d3c9b50f19704552f23b4eaea1fc082fdd82c63429a6506446cbd8737823da3"}, + {file = "aiohttp-3.9.3-cp311-cp311-musllinux_1_1_ppc64le.whl", hash = "sha256:f033d80bc6283092613882dfe40419c6a6a1527e04fc69350e87a9df02bbc283"}, + {file = "aiohttp-3.9.3-cp311-cp311-musllinux_1_1_s390x.whl", hash = "sha256:2c895a656dd7e061b2fd6bb77d971cc38f2afc277229ce7dd3552de8313a483e"}, + {file = "aiohttp-3.9.3-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:1f5a71d25cd8106eab05f8704cd9167b6e5187bcdf8f090a66c6d88b634802b4"}, + {file = "aiohttp-3.9.3-cp311-cp311-win32.whl", hash = "sha256:50fca156d718f8ced687a373f9e140c1bb765ca16e3d6f4fe116e3df7c05b2c5"}, + {file = "aiohttp-3.9.3-cp311-cp311-win_amd64.whl", hash = "sha256:5fe9ce6c09668063b8447f85d43b8d1c4e5d3d7e92c63173e6180b2ac5d46dd8"}, + {file = "aiohttp-3.9.3-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:38a19bc3b686ad55804ae931012f78f7a534cce165d089a2059f658f6c91fa60"}, + {file = "aiohttp-3.9.3-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:770d015888c2a598b377bd2f663adfd947d78c0124cfe7b959e1ef39f5b13869"}, + {file = "aiohttp-3.9.3-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:ee43080e75fc92bf36219926c8e6de497f9b247301bbf88c5c7593d931426679"}, + {file = "aiohttp-3.9.3-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:52df73f14ed99cee84865b95a3d9e044f226320a87af208f068ecc33e0c35b96"}, + {file = "aiohttp-3.9.3-cp312-cp312-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:dc9b311743a78043b26ffaeeb9715dc360335e5517832f5a8e339f8a43581e4d"}, + {file = "aiohttp-3.9.3-cp312-cp312-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:b955ed993491f1a5da7f92e98d5dad3c1e14dc175f74517c4e610b1f2456fb11"}, + {file = "aiohttp-3.9.3-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:504b6981675ace64c28bf4a05a508af5cde526e36492c98916127f5a02354d53"}, + {file = "aiohttp-3.9.3-cp312-cp312-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:a6fe5571784af92b6bc2fda8d1925cccdf24642d49546d3144948a6a1ed58ca5"}, + {file = "aiohttp-3.9.3-cp312-cp312-musllinux_1_1_aarch64.whl", hash = "sha256:ba39e9c8627edc56544c8628cc180d88605df3892beeb2b94c9bc857774848ca"}, + {file = "aiohttp-3.9.3-cp312-cp312-musllinux_1_1_i686.whl", hash = "sha256:e5e46b578c0e9db71d04c4b506a2121c0cb371dd89af17a0586ff6769d4c58c1"}, + {file = "aiohttp-3.9.3-cp312-cp312-musllinux_1_1_ppc64le.whl", hash = "sha256:938a9653e1e0c592053f815f7028e41a3062e902095e5a7dc84617c87267ebd5"}, + {file = "aiohttp-3.9.3-cp312-cp312-musllinux_1_1_s390x.whl", hash = "sha256:c3452ea726c76e92f3b9fae4b34a151981a9ec0a4847a627c43d71a15ac32aa6"}, + {file = "aiohttp-3.9.3-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:ff30218887e62209942f91ac1be902cc80cddb86bf00fbc6783b7a43b2bea26f"}, + {file = "aiohttp-3.9.3-cp312-cp312-win32.whl", hash = "sha256:38f307b41e0bea3294a9a2a87833191e4bcf89bb0365e83a8be3a58b31fb7f38"}, + {file = "aiohttp-3.9.3-cp312-cp312-win_amd64.whl", hash = "sha256:b791a3143681a520c0a17e26ae7465f1b6f99461a28019d1a2f425236e6eedb5"}, + {file = "aiohttp-3.9.3-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:0ed621426d961df79aa3b963ac7af0d40392956ffa9be022024cd16297b30c8c"}, + {file = "aiohttp-3.9.3-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:7f46acd6a194287b7e41e87957bfe2ad1ad88318d447caf5b090012f2c5bb528"}, + {file = "aiohttp-3.9.3-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:feeb18a801aacb098220e2c3eea59a512362eb408d4afd0c242044c33ad6d542"}, + {file = "aiohttp-3.9.3-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f734e38fd8666f53da904c52a23ce517f1b07722118d750405af7e4123933511"}, + {file = "aiohttp-3.9.3-cp38-cp38-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:b40670ec7e2156d8e57f70aec34a7216407848dfe6c693ef131ddf6e76feb672"}, + {file = "aiohttp-3.9.3-cp38-cp38-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:fdd215b7b7fd4a53994f238d0f46b7ba4ac4c0adb12452beee724ddd0743ae5d"}, + {file = "aiohttp-3.9.3-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:017a21b0df49039c8f46ca0971b3a7fdc1f56741ab1240cb90ca408049766168"}, + {file = "aiohttp-3.9.3-cp38-cp38-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:e99abf0bba688259a496f966211c49a514e65afa9b3073a1fcee08856e04425b"}, + {file = "aiohttp-3.9.3-cp38-cp38-musllinux_1_1_aarch64.whl", hash = "sha256:648056db9a9fa565d3fa851880f99f45e3f9a771dd3ff3bb0c048ea83fb28194"}, + {file = "aiohttp-3.9.3-cp38-cp38-musllinux_1_1_i686.whl", hash = "sha256:8aacb477dc26797ee089721536a292a664846489c49d3ef9725f992449eda5a8"}, + {file = "aiohttp-3.9.3-cp38-cp38-musllinux_1_1_ppc64le.whl", hash = "sha256:522a11c934ea660ff8953eda090dcd2154d367dec1ae3c540aff9f8a5c109ab4"}, + {file = "aiohttp-3.9.3-cp38-cp38-musllinux_1_1_s390x.whl", hash = "sha256:5bce0dc147ca85caa5d33debc4f4d65e8e8b5c97c7f9f660f215fa74fc49a321"}, + {file = "aiohttp-3.9.3-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:4b4af9f25b49a7be47c0972139e59ec0e8285c371049df1a63b6ca81fdd216a2"}, + {file = "aiohttp-3.9.3-cp38-cp38-win32.whl", hash = "sha256:298abd678033b8571995650ccee753d9458dfa0377be4dba91e4491da3f2be63"}, + {file = "aiohttp-3.9.3-cp38-cp38-win_amd64.whl", hash = "sha256:69361bfdca5468c0488d7017b9b1e5ce769d40b46a9f4a2eed26b78619e9396c"}, + {file = "aiohttp-3.9.3-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:0fa43c32d1643f518491d9d3a730f85f5bbaedcbd7fbcae27435bb8b7a061b29"}, + {file = "aiohttp-3.9.3-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:835a55b7ca49468aaaac0b217092dfdff370e6c215c9224c52f30daaa735c1c1"}, + {file = "aiohttp-3.9.3-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:06a9b2c8837d9a94fae16c6223acc14b4dfdff216ab9b7202e07a9a09541168f"}, + {file = "aiohttp-3.9.3-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:abf151955990d23f84205286938796c55ff11bbfb4ccfada8c9c83ae6b3c89a3"}, + {file = "aiohttp-3.9.3-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:59c26c95975f26e662ca78fdf543d4eeaef70e533a672b4113dd888bd2423caa"}, + {file = "aiohttp-3.9.3-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:f95511dd5d0e05fd9728bac4096319f80615aaef4acbecb35a990afebe953b0e"}, + {file = "aiohttp-3.9.3-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:595f105710293e76b9dc09f52e0dd896bd064a79346234b521f6b968ffdd8e58"}, + {file = "aiohttp-3.9.3-cp39-cp39-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:c7c8b816c2b5af5c8a436df44ca08258fc1a13b449393a91484225fcb7545533"}, + {file = "aiohttp-3.9.3-cp39-cp39-musllinux_1_1_aarch64.whl", hash = "sha256:f1088fa100bf46e7b398ffd9904f4808a0612e1d966b4aa43baa535d1b6341eb"}, + {file = "aiohttp-3.9.3-cp39-cp39-musllinux_1_1_i686.whl", hash = "sha256:f59dfe57bb1ec82ac0698ebfcdb7bcd0e99c255bd637ff613760d5f33e7c81b3"}, + {file = "aiohttp-3.9.3-cp39-cp39-musllinux_1_1_ppc64le.whl", hash = "sha256:361a1026c9dd4aba0109e4040e2aecf9884f5cfe1b1b1bd3d09419c205e2e53d"}, + {file = "aiohttp-3.9.3-cp39-cp39-musllinux_1_1_s390x.whl", hash = "sha256:363afe77cfcbe3a36353d8ea133e904b108feea505aa4792dad6585a8192c55a"}, + {file = "aiohttp-3.9.3-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:8e2c45c208c62e955e8256949eb225bd8b66a4c9b6865729a786f2aa79b72e9d"}, + {file = "aiohttp-3.9.3-cp39-cp39-win32.whl", hash = "sha256:f7217af2e14da0856e082e96ff637f14ae45c10a5714b63c77f26d8884cf1051"}, + {file = "aiohttp-3.9.3-cp39-cp39-win_amd64.whl", hash = "sha256:27468897f628c627230dba07ec65dc8d0db566923c48f29e084ce382119802bc"}, + {file = "aiohttp-3.9.3.tar.gz", hash = "sha256:90842933e5d1ff760fae6caca4b2b3edba53ba8f4b71e95dacf2818a2aca06f7"}, ] [package.dependencies] @@ -128,9 +128,6 @@ files = [ {file = "aioitertools-0.11.0.tar.gz", hash = "sha256:42c68b8dd3a69c2bf7f2233bf7df4bb58b557bca5252ac02ed5187bbc67d6831"}, ] -[package.dependencies] -typing_extensions = {version = ">=4.0", markers = "python_version < \"3.10\""} - [[package]] name = "aiosignal" version = "1.3.1" @@ -359,10 +356,7 @@ files = [ [package.dependencies] jmespath = ">=0.7.1,<2.0.0" python-dateutil = ">=2.1,<3.0.0" -urllib3 = [ - {version = ">=1.25.4,<1.27", markers = "python_version < \"3.10\""}, - {version = ">=1.25.4,<2.1", markers = "python_version >= \"3.10\""}, -] +urllib3 = {version = ">=1.25.4,<2.1", markers = "python_version >= \"3.10\""} [package.extras] crt = ["awscrt (==0.19.19)"] @@ -832,6 +826,43 @@ files = [ [package.dependencies] pyyaml = "*" +[[package]] +name = "ducc0" +version = "0.31.0" +description = "Distinctly useful code collection: contains efficient algorithms for Fast Fourier (and related) transforms, spherical harmonic transforms involving very general spherical grids, gridding/degridding tools for radio interferometry, 4pi spherical convolution operators and much more." +optional = false +python-versions = ">=3.7" +files = [ + {file = "ducc0-0.31.0-cp310-cp310-macosx_10_14_x86_64.whl", hash = "sha256:268c70d31047afaa07ec3db336d0b0de148bdaacf7a7f2261f99cfcda1572108"}, + {file = "ducc0-0.31.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:e5e3d2be28693d18ea8d1ea8fe55ad8e526ab037ce472c1fee4772ee738d2c27"}, + {file = "ducc0-0.31.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:a07dc22c3e36cd1f5b1090e57825ec9de24785482e05d38e49171664743c67e1"}, + {file = "ducc0-0.31.0-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:65dd2cff8cbf55eedd65b9bab500bc63da37a28753697d9c296f3f9a377b5754"}, + {file = "ducc0-0.31.0-cp310-cp310-win_amd64.whl", hash = "sha256:795e8f8bf131b7b62b97c30eb0e7dd523577e16964b66b7fb7cad5dc02603711"}, + {file = "ducc0-0.31.0-cp311-cp311-macosx_10_14_x86_64.whl", hash = "sha256:35ee18e8a184692b1a4a0b853cbecff2884d97e3951a35e232395b72fd3c4d6f"}, + {file = "ducc0-0.31.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:7fed867befacd39866cf4d050cee71a4fca3aa59084c02234913d7f54b43c196"}, + {file = "ducc0-0.31.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e270f62cdd0184b85a5f4084f160ee6315ffc62ba5eb02ce59344d62db7748ca"}, + {file = "ducc0-0.31.0-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:7b6ac08bde0c3c303b190a8847bffafad47e328af0df66ae286ea4435ba6e2b5"}, + {file = "ducc0-0.31.0-cp311-cp311-win_amd64.whl", hash = "sha256:2649faff4479746ad0bdd42b061fc2992314542239a22bb597a5bb443bdd2c28"}, + {file = "ducc0-0.31.0-cp37-cp37m-macosx_10_14_x86_64.whl", hash = "sha256:16ea97198ab37262d0e2ede1956b0bfe58d15a37c278fc280de4e362ffdd6663"}, + {file = "ducc0-0.31.0-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:01b18372840c3650561a461d906fb568bc944639f59db08d32846d5edb2d5346"}, + {file = "ducc0-0.31.0-cp37-cp37m-musllinux_1_1_x86_64.whl", hash = "sha256:b4da58fb5eff415dc1c3e87e7272846abff935707536a8efeaf14399475c4c8f"}, + {file = "ducc0-0.31.0-cp37-cp37m-win_amd64.whl", hash = "sha256:7c34ac70c2b47de312dc3eb55af2208ea63dbc174a3c7e66587e3ffb67c1b98e"}, + {file = "ducc0-0.31.0-cp38-cp38-macosx_10_14_x86_64.whl", hash = "sha256:4e86ac3251bdd89e65c24c55f26b4e363f1502488d33e128f4e35c363ded81ff"}, + {file = "ducc0-0.31.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:578b36755cf1b8b37bd2b1e6ffc28a73846eaed3235309cc35fef2275fc37c9a"}, + {file = "ducc0-0.31.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:2d60392072deef7ea0b1e6f79bd321c8f7988b132c39b1a6d137a9b3388faaa6"}, + {file = "ducc0-0.31.0-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:c8f742d934e5afcde40835c759b901b8b50b659ea194d766299ae83be3adedda"}, + {file = "ducc0-0.31.0-cp38-cp38-win_amd64.whl", hash = "sha256:2ba8e5047e8d6420ab6666a43b32d617690f84579872e88e79297ce1785c3dcf"}, + {file = "ducc0-0.31.0-cp39-cp39-macosx_10_14_x86_64.whl", hash = "sha256:513ba3c5a8d27d24717eff1a549f5fca5c75a94d9218c572a15efa59104f20ee"}, + {file = "ducc0-0.31.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:a15de0768967b7e80dc1ed227b69f089d9c3a7edd7523b83049965107d1ea508"}, + {file = "ducc0-0.31.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:c821334dde973acc818ca6f0ca0ef239812e237f6f738842795fe7afaf71f96d"}, + {file = "ducc0-0.31.0-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:d4e6432bce517ba8f92d4e085edfc6868228541acf117bb0ee495549492b7c3f"}, + {file = "ducc0-0.31.0-cp39-cp39-win_amd64.whl", hash = "sha256:81144c963576978cdcc783f922c421538b918f75ba59d7c312c9ffb296c3b345"}, + {file = "ducc0-0.31.0.tar.gz", hash = "sha256:ff2387f292d33de9fc8804df6a957f50a11474ceef65532d37afa0a4c333c9e5"}, +] + +[package.dependencies] +numpy = ">=1.17.0" + [[package]] name = "exceptiongroup" version = "1.2.0" @@ -1083,24 +1114,6 @@ docs = ["furo", "jaraco.packaging (>=9.3)", "jaraco.tidelift (>=1.4)", "rst.link perf = ["ipython"] testing = ["flufl.flake8", "importlib-resources (>=1.3)", "packaging", "pyfakefs", "pytest (>=6)", "pytest-black (>=0.3.7)", "pytest-checkdocs (>=2.4)", "pytest-cov", "pytest-enabler (>=2.2)", "pytest-mypy (>=0.9.1)", "pytest-perf (>=0.9.2)", "pytest-ruff"] -[[package]] -name = "importlib-resources" -version = "6.1.1" -description = "Read resources from Python packages" -optional = false -python-versions = ">=3.8" -files = [ - {file = "importlib_resources-6.1.1-py3-none-any.whl", hash = "sha256:e8bf90d8213b486f428c9c39714b920041cb02c184686a3dee24905aaa8105d6"}, - {file = "importlib_resources-6.1.1.tar.gz", hash = "sha256:3893a00122eafde6894c59914446a512f728a0c1a45f9bb9b63721b6bacf0b4a"}, -] - -[package.dependencies] -zipp = {version = ">=3.1.0", markers = "python_version < \"3.10\""} - -[package.extras] -docs = ["furo", "jaraco.packaging (>=9.3)", "jaraco.tidelift (>=1.4)", "rst.linker (>=1.9)", "sphinx (<7.2.5)", "sphinx (>=3.5)", "sphinx-lint"] -testing = ["pytest (>=6)", "pytest-black (>=0.3.7)", "pytest-checkdocs (>=2.4)", "pytest-cov", "pytest-enabler (>=2.2)", "pytest-mypy (>=0.9.1)", "pytest-ruff", "zipp (>=3.17)"] - [[package]] name = "iniconfig" version = "2.0.0" @@ -1425,7 +1438,6 @@ files = [ contourpy = ">=1.0.1" cycler = ">=0.10" fonttools = ">=4.22.0" -importlib-resources = {version = ">=3.2.0", markers = "python_version < \"3.10\""} kiwisolver = ">=1.3.1" numpy = ">=1.21,<2" packaging = ">=20.0" @@ -1433,6 +1445,23 @@ pillow = ">=8" pyparsing = ">=2.3.1" python-dateutil = ">=2.7" +[[package]] +name = "mpmath" +version = "1.3.0" +description = "Python library for arbitrary-precision floating-point arithmetic" +optional = false +python-versions = "*" +files = [ + {file = "mpmath-1.3.0-py3-none-any.whl", hash = "sha256:a0b2b9fe80bbcd81a6647ff13108738cfb482d481d826cc0e02f5b35e5c88d2c"}, + {file = "mpmath-1.3.0.tar.gz", hash = "sha256:7a28eb2a9774d00c7bc92411c19a89209d5da7c4c9a9e227be8330a23a25b91f"}, +] + +[package.extras] +develop = ["codecov", "pycodestyle", "pytest (>=4.6)", "pytest-cov", "wheel"] +docs = ["sphinx"] +gmpy = ["gmpy2 (>=2.1.0a4)"] +tests = ["pytest (>=4.6)"] + [[package]] name = "msgpack" version = "1.0.7" @@ -1599,6 +1628,32 @@ six = "*" testing = ["astroid (>=1.5.3,<1.6.0)", "astroid (>=2.0)", "coverage", "pylint (>=1.7.2,<1.8.0)", "pylint (>=2.3.1,<2.4.0)", "pytest"] yaml = ["PyYAML (>=5.1.0)"] +[[package]] +name = "nifty8" +version = "8.4" +description = "Library for signal inference algorithms that operate regardless of the underlying grids and their resolutions." +optional = false +python-versions = ">=3.10" +files = [] +develop = false + +[package.dependencies] +numpy = ">=1.23" +scipy = ">=1.9.0" + +[package.extras] +doc = ["jupyter", "jupytext", "pydata-sphinx-theme", "sphinx"] +full = ["astropy", "ducc0 (>=0.27.0)", "healpy", "jax", "jaxlib", "jupyter", "jupytext", "pydata-sphinx-theme", "sphinx"] +native = ["ducc0 (>=0.27.0)"] +re = ["jax", "jaxlib"] +util = ["astropy", "healpy"] + +[package.source] +type = "git" +url = "https://gitlab.mpcdf.mpg.de/ift/nifty.git" +reference = "NIFTy_8" +resolved_reference = "dda7464c34bf28ea08f0e30763a3d434ebdae499" + [[package]] name = "numba" version = "0.58.1" @@ -2432,6 +2487,20 @@ pydantic = ">=1.10.2,<2.0.0" pyparsing = ">=3.0.9,<4.0.0" rich = ">=12.6.0,<13.0.0" +[[package]] +name = "sympy" +version = "1.12" +description = "Computer algebra system (CAS) in Python" +optional = false +python-versions = ">=3.8" +files = [ + {file = "sympy-1.12-py3-none-any.whl", hash = "sha256:c3588cd4295d0c0f603d0f2ae780587e64e2efeedb3521e46b9bb1d08d184fa5"}, + {file = "sympy-1.12.tar.gz", hash = "sha256:ebf595c8dac3e0fdc4152c51878b498396ec7f30e7a914d6071e674d49420fb8"}, +] + +[package.dependencies] +mpmath = ">=0.19" + [[package]] name = "tabulate" version = "0.8.10" @@ -2560,22 +2629,6 @@ files = [ {file = "Unidecode-1.3.8.tar.gz", hash = "sha256:cfdb349d46ed3873ece4586b96aa75258726e2fa8ec21d6f00a591d98806c2f4"}, ] -[[package]] -name = "urllib3" -version = "1.26.18" -description = "HTTP library with thread-safe connection pooling, file post, and more." -optional = false -python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*, !=3.5.*" -files = [ - {file = "urllib3-1.26.18-py2.py3-none-any.whl", hash = "sha256:34b97092d7e0a3a8cf7cd10e386f401b3737364026c45e622aa02903dffe0f07"}, - {file = "urllib3-1.26.18.tar.gz", hash = "sha256:f8ecc1bba5667413457c529ab955bf8c67b45db799d159066261719e328580a0"}, -] - -[package.extras] -brotli = ["brotli (==1.0.9)", "brotli (>=1.0.9)", "brotlicffi (>=0.8.0)", "brotlipy (>=0.6.0)"] -secure = ["certifi", "cryptography (>=1.3.4)", "idna (>=2.0.0)", "ipaddress", "pyOpenSSL (>=0.14)", "urllib3-secure-extra"] -socks = ["PySocks (>=1.5.6,!=1.5.7,<2.0)"] - [[package]] name = "urllib3" version = "2.0.7" @@ -2881,7 +2934,10 @@ files = [ docs = ["furo", "jaraco.packaging (>=9.3)", "jaraco.tidelift (>=1.4)", "rst.linker (>=1.9)", "sphinx (<7.2.5)", "sphinx (>=3.5)", "sphinx-lint"] testing = ["big-O", "jaraco.functools", "jaraco.itertools", "more-itertools", "pytest (>=6)", "pytest-black (>=0.3.7)", "pytest-checkdocs (>=2.4)", "pytest-cov", "pytest-enabler (>=2.2)", "pytest-ignore-flaky", "pytest-mypy (>=0.9.1)", "pytest-ruff"] +[extras] +degrid = ["ducc0", "sympy"] + [metadata] lock-version = "2.0" -python-versions = "^3.9" -content-hash = "59eb838a688ad3f3c25ab91cafc12d15fa6d098d3b339458bcabc294531a009e" +python-versions = "^3.10" +content-hash = "e4be9fd30fc0d021e9bf9eea45be7a46d96de81351454c29fb867da8a47e6c59" diff --git a/pyproject.toml b/pyproject.toml index 5a350c6b..ab0a7b3e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,6 +44,8 @@ ducc0 = "^0.31.0" sympy = "^1.12" nifty8 = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch = "NIFTy_8"} +[tool.poetry.extras] +degrid = ["ducc0", "sympy"] [tool.poetry.scripts] goquartical = 'quartical.executor:execute' From dc56fb55d65c4313175ceb09a2bc735c6b5ea9f4 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 30 Jan 2024 14:42:16 +0200 Subject: [PATCH 62/63] xds_from_fragment_ms in ms_handler --- quartical/data_handling/ms_handler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py index 57cdce45..d6051aba 100644 --- a/quartical/data_handling/ms_handler.py +++ b/quartical/data_handling/ms_handler.py @@ -97,7 +97,7 @@ def read_xds_list(model_columns, ms_opts): if ms_opts.data_column not in known_data_cols: schema[ms_opts.data_column] = {'dims': ('chan', 'corr')} - data_xds_list = xds_from_storage_ms( + data_xds_list = xds_from_ms_fragment( ms_opts.path, columns=columns, index_cols=("TIME",), From 829948f13e8404422492a8c750bcf7f5db4869ab Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 30 Jan 2024 15:24:19 +0200 Subject: [PATCH 63/63] fix unmerged conflicts --- quartical/apps/plotter.py | 20 -------------------- quartical/apps/summary.py | 4 ---- 2 files changed, 24 deletions(-) diff --git a/quartical/apps/plotter.py b/quartical/apps/plotter.py index b3e6785b..35f18977 100644 --- a/quartical/apps/plotter.py +++ b/quartical/apps/plotter.py @@ -216,16 +216,9 @@ def _plot(group, xds, args): # NOTE: This mututates the data variables in place. data = xds[args.plot_var].values -<<<<<<< HEAD - # flags = xds[args.flag_var].values - # data[np.where(flags)] = np.nan # Set flagged values to nan (not plotted). - # xds['flagged_data'] = np.where(flags, data, np.nan) - # xds = xds.drop_vars(args.flag_var) # No more use for flags. -======= flags = xds[args.flag_var].values data[np.where(flags)] = np.nan # Set flagged values to nan (not plotted). xds = xds.drop_vars(args.flag_var) # No more use for flags. ->>>>>>> main # Construct list of lists containing axes over which we iterate i.e. # produce a plot per combination of these values. @@ -251,10 +244,6 @@ def _plot(group, xds, args): fig, ax = plt.subplots(figsize=args.fig_size) -<<<<<<< HEAD - # import ipdb; ipdb.set_trace() -======= ->>>>>>> main for ia in product(*iter_axes_itr): sel = {ax: val for ax, val in zip(args.iter_axes, ia)} @@ -311,10 +300,6 @@ def _plot(group, xds, args): fig.savefig( f"{args.output_path.full_path}/{subdir_path}/{fig_name}.png", -<<<<<<< HEAD - dpi=250, -======= ->>>>>>> main bbox_inches="tight" # SLOW, but slightly unavoidable. ) @@ -335,12 +320,7 @@ def plot(): xdsl = xds_from_zarr(gain_path) # Select only the necessary fields for plotting on each dataset. -<<<<<<< HEAD - # xdsl = [xds[[args.plot_var, args.flag_var]] for xds in xdsl] - xdsl = [xds[[args.plot_var]] for xds in xdsl] -======= xdsl = [xds[[args.plot_var, args.flag_var]] for xds in xdsl] ->>>>>>> main # Partitioned dictionary of xarray.Datasets. xdsd = to_plot_dict(xdsl, args.iter_attrs) diff --git a/quartical/apps/summary.py b/quartical/apps/summary.py index b6abfb9d..716b0653 100644 --- a/quartical/apps/summary.py +++ b/quartical/apps/summary.py @@ -119,11 +119,7 @@ def field_info(path): field_xds = xds_from_table_fragment(path + "::FIELD")[0] -<<<<<<< HEAD - field_ids = list(range(field_xds.dims['row'])) -======= field_ids = list(range(field_xds.sizes['row'])) ->>>>>>> main source_ids = [i for i in field_xds.SOURCE_ID.values] names = [n for n in field_xds.NAME.values] phase_dirs = [pd for pd in field_xds.PHASE_DIR.values]