From fdad13812bc1323bbb03a5dc4120e850f406ac37 Mon Sep 17 00:00:00 2001 From: Emanuele Zanoni Date: Tue, 13 Jan 2026 18:42:58 +0100 Subject: [PATCH 1/6] fix(bayesian_listener): apply Ruff fixes (except ERA and D warnings) --- bayesian_listener/bayesian_listener.py | 205 ++++++++++++++++++------- 1 file changed, 147 insertions(+), 58 deletions(-) diff --git a/bayesian_listener/bayesian_listener.py b/bayesian_listener/bayesian_listener.py index 2ba29d0..95fe4a0 100644 --- a/bayesian_listener/bayesian_listener.py +++ b/bayesian_listener/bayesian_listener.py @@ -26,11 +26,13 @@ def __init__(self, sofa): self.fs = int(self.sofa_data.Data_SamplingRate) self.coords = Coordinates(sofa_file = sofa) # noise parameters - self.parameters = dict(sigma_itd = .569, - sigma_ild = 0.75, - sigma_spectral = 4, - sigma_prior = 11.5, - sigma_motor = 12) + self.parameters = { + "sigma_itd": 0.569, + "sigma_ild": 0.75, + "sigma_spectral": 4, + "sigma_prior": 11.5, + "sigma_motor": 12, + } @property def parameters(self): @@ -41,29 +43,41 @@ def parameters(self, value): if not isinstance(value, dict): raise ValueError("Parameters must be a dictionary.") # check if all parameters are present - for key in ['sigma_itd', 'sigma_ild', 'sigma_spectral', 'sigma_prior', 'sigma_motor']: + for key in [ + 'sigma_itd', + 'sigma_ild', + 'sigma_spectral', + 'sigma_prior', + 'sigma_motor', + ]: if key not in value: raise ValueError(f"Missing parameter: {key}") self._parameters = value - # interpolate to a "uniform" spherical grid coz usually HRTF do not have a uniform grid + # interpolate to a "uniform" spherical grid + # coz usually HRTF do not have a uniform grid def interpolate(self, interpolation='SH'): - model = BayesianListener.__new__(BayesianListener) # Create empty instance + # Create empty instance + model = BayesianListener.__new__(BayesianListener) # Resample all cues in a single call cues_list = [ self.itd, self.ild, self.spectral_cues[:, :, 0], - self.spectral_cues[:, :, 1] + self.spectral_cues[:, :, 1], ] - resampled_cues, coords_new = resample.resample(cues_list, self.coords, method=interpolation) + resampled_cues, coords_new = resample.resample(cues_list, + self.coords, + method=interpolation) # Unpack results model.itd = resampled_cues[0] model.ild = resampled_cues[1] - model.spectral_cues = np.stack([resampled_cues[2], resampled_cues[3]], axis=-1) + model.spectral_cues = np.stack([resampled_cues[2], + resampled_cues[3]], + axis=-1) model.coords = coords_new model.freqs = self.freqs # model.coords.plot(model.spectral_cues[:, 5, 0]) @@ -72,17 +86,28 @@ def interpolate(self, interpolation='SH'): return model # prepare - def prepare_features(self, spectral_range=[7e2, 18e3], interpolation='SH', use_cache=True, force_recompute=False): - """delegates to cached or direct computation of spatial features and templates""" + def prepare_features(self, + spectral_range=[7e2, 18e3], + interpolation='SH', + use_cache=True, + force_recompute=False): + """ + delegates to cached or direct computation of spatial features and + templates. + """ assert(self.sofa_file is not None) if use_cache: - return self._load_or_compute_features(spectral_range, interpolation, force_recompute) + return self._load_or_compute_features(spectral_range, + interpolation, + force_recompute) else: return self._compute_features(spectral_range, interpolation) # pre-compute features for faster inference - def _compute_features(self, spectral_range = [7e2, 18e3], interpolation='SH'): + def _compute_features(self, + spectral_range = [7e2, 18e3], + interpolation='SH'): # normalize hrirs to frontal position _, idx = self.coords.find(Coordinates(positions=np.array([1, 0, 0]))) hrirs_temp = self.hrir / np.max(np.abs(self.hrir[idx])) @@ -92,12 +117,14 @@ def _compute_features(self, spectral_range = [7e2, 18e3], interpolation='SH'): # ITD itd = utils.itdestimator(hrirs_temp, fs=self.fs) - self.itd = np.sign(itd) * ((np.log(a + b * np.abs(itd)) - np.log(a)) / b) + self.itd = np.sign(itd) * ((np.log(a + b*np.abs(itd)) - np.log(a)) / b) # ILD self.ild = np.ones_like(self.itd) - self.ild[:, 0] = (utils.mag2db(np.sqrt(np.mean(hrirs_temp[:, 0, :]**2, axis=1))) - - utils.mag2db(np.sqrt(np.mean(hrirs_temp[:, 1, :]**2, axis=1)))) + self.ild[:, 0] = ( + utils.mag2db(np.sqrt(np.mean(hrirs_temp[:, 0, :]**2,axis=1))) - + utils.mag2db(np.sqrt(np.mean(hrirs_temp[:, 1, :]**2, axis=1))) + ) # compute spatial features # -------- padding to account for longer filter responses -------- @@ -109,7 +136,8 @@ def _compute_features(self, spectral_range = [7e2, 18e3], interpolation='SH'): if time_len < target_samples: pad_samples = target_samples - time_len - pad_mat = np.zeros((dir_len, ear_len, pad_samples), dtype=hrirs_temp.dtype) + pad_mat = np.zeros((dir_len, ear_len, pad_samples), + dtype=hrirs_temp.dtype) hrirs_temp = np.concatenate([hrirs_temp, pad_mat], axis=2) # generate gammatone filterbank @@ -117,7 +145,8 @@ def _compute_features(self, spectral_range = [7e2, 18e3], interpolation='SH'): B, A, *_ = utils.gammatone(self.freqs, fs=self.fs) # Preallocate output array (float, since we take 2*real(...)) - hrirs_filt = np.zeros((len(self.freqs), *hrirs_temp.shape), dtype=float) + hrirs_filt = np.zeros((len(self.freqs), *hrirs_temp.shape), + dtype=float) # Parallel gammatone filtering def apply_filter(i): @@ -141,17 +170,23 @@ def apply_filter(i): hrirs_filt = np.sqrt(np.maximum(hrirs_filt, 0)) # average over time -> spectral amplitude - rms = np.sqrt(np.mean(hrirs_filt**2, axis=-1)) # (n_freqs, n_dirs, n_ears) - spectral_amps = utils.mag2db(rms).transpose(1, 0, 2) # -> (n_dirs, n_freqs, n_ears) + # (n_freqs, n_dirs, n_ears) + rms = np.sqrt(np.mean(hrirs_filt**2, axis=-1)) + # -> (n_dirs, n_freqs, n_ears) + spectral_amps = utils.mag2db(rms).transpose(1, 0, 2) self.spectral_cues = spectral_amps # prepare templates on uniform grid self.template = self.interpolate(interpolation) - def _load_or_compute_features(self, spectral_range=[7e2, 18e3], interpolation='SH', force_recompute=False): + def _load_or_compute_features(self, + spectral_range=[7e2, 18e3], + interpolation='SH', + force_recompute=False): """ - Preprocess with caching: load from cache if available, otherwise compute and save. + Preprocess with caching: load from cache if available, + otherwise compute and save. Parameters ---------- @@ -174,12 +209,15 @@ def _load_or_compute_features(self, spectral_range=[7e2, 18e3], interpolation='S # Define what attributes to cache/restore cache_attributes = [ 'itd', 'ild', 'freqs', 'spectral_cues', - 'coords', 'parameters', 'template' + 'coords', 'parameters', 'template', ] # ========== Try to load from cache ========== if not force_recompute and self.sofa_file is not None: - cached_data = utils.load_from_cache(cache_dir, self.sofa_file, cache_attributes, interpolation) + cached_data = utils.load_from_cache(cache_dir, + self.sofa_file, + cache_attributes, + interpolation) if cached_data is not None: # Restore cached attributes @@ -191,23 +229,38 @@ def _load_or_compute_features(self, spectral_range=[7e2, 18e3], interpolation='S # ========== Compute features ========== print("→ Computing features...") - self._compute_features(spectral_range=spectral_range, interpolation=interpolation) + self._compute_features(spectral_range=spectral_range, + interpolation=interpolation) print("✓ Feature preparation complete") # ========== Save to cache ========== if self.sofa_file is not None: # Prepare data to cache - cache_data = {attr: getattr(self, attr) for attr in cache_attributes} - utils.save_to_cache(cache_dir, self.sofa_file, cache_data, interpolation) + cache_data = { + attr: getattr(self, attr) for attr in cache_attributes + } + utils.save_to_cache(cache_dir, + self.sofa_file, + cache_data, + interpolation) # return internal representation def represent(self): - bcue = np.hstack([self.itd, self.ild]) - scue = np.hstack([self.spectral_cues[:, :, 0], self.spectral_cues[:, :, 1]]) - return np.hstack([bcue, scue]) + bcue = np.hstack([self.itd, + self.ild]) + + scue = np.hstack([self.spectral_cues[:, :, 0], + self.spectral_cues[:, :, 1]]) + + return np.hstack([bcue, + scue]) - def infer(self, target = None, repetitions = 50, seed = None, prior = 'horizontal'): - np.random.seed(seed) + def infer(self, + target = None, + repetitions = 50, + seed = None, + prior = 'horizontal'): + rng = np.random.default_rng(seed) # prepare features # use original HRIR if no target is provided @@ -224,8 +277,11 @@ def infer(self, target = None, repetitions = 50, seed = None, prior = 'horizonta template_feat = self.template.represent() sigmas = self.parameters - sigma = np.block(np.diag(np.hstack([sigmas["sigma_itd"]**2, sigmas["sigma_ild"]**2, - np.repeat(sigmas["sigma_spectral"]**2, self.freqs.shape[0]*2)]))) + sigma = np.block(np.diag(np.hstack( + [sigmas["sigma_itd"]**2, + sigmas["sigma_ild"]**2, + np.repeat(sigmas["sigma_spectral"]**2, self.freqs.shape[0]*2), + ]))) # the following code is needed to speed up multiple_logpdfs_vec_input # since here the covariance matrix is constant NumPy broadcasts `eigh`. @@ -233,7 +289,8 @@ def infer(self, target = None, repetitions = 50, seed = None, prior = 'horizonta # Compute the log determinants across the second axis. logdet = np.sum(np.log(vals)) - # Invert the eigenvalues and add a dimension to `valsinvs` so that NumPy broadcasts appropriately. + # Invert the eigenvalues and add a dimension to `valsinvs` + # so that NumPy broadcasts appropriately. Us = vecs * np.sqrt(1./vals)[:, None] # Prior computation @@ -242,20 +299,28 @@ def infer(self, target = None, repetitions = 50, seed = None, prior = 'horizonta # Uniform prior: all directions equally likely prior = np.ones(template_feat.shape[0]) elif prior == 'horizontal': - # Horizontal bias prior: Gaussian centered on horizontal plane (elevation = 0°) + # Horizontal bias prior: + # Gaussian centered on horizontal plane (elevation = 0°) sph = self.template.coords.convert('spherical') - prior = np.exp(-0.5 * (np.rad2deg(sph[:, 1]) / sigmas["sigma_prior"])**2) + prior = np.exp( + -0.5 * (np.rad2deg(sph[:, 1]) / sigmas["sigma_prior"])**2, + ) else: - raise ValueError(f"Unknown prior: {prior}. Use 'uniform', 'horizontal', or numpy array") + raise ValueError( + f"Unknown prior: {prior}. " + f"Use 'uniform', 'horizontal', or numpy array") # Normalize to sum to 1 (valid probability distribution) prior /= np.sum(prior) elif isinstance(prior, np.ndarray): # Custom prior provided as array if prior.shape[0] != template_feat.shape[0]: - raise ValueError(f"Prior shape mismatch: {prior.shape[0]} vs {template_feat.shape[0]}") + raise ValueError( + f"Prior shape mismatch: " + f"{prior.shape[0]} vs {template_feat.shape[0]}") prior = prior / np.sum(prior) # normalize else: - raise TypeError("Prior must be str ('uniform', 'horizontal') or numpy array") + raise TypeError( + "Prior must be str ('uniform', 'horizontal') or numpy array") # self.template.coords.plot(prior) @@ -267,44 +332,56 @@ def infer(self, target = None, repetitions = 50, seed = None, prior = 'horizonta L = np.linalg.cholesky(sigma) # L @ L.T = sigma for t in range(target_num): ts = np.tile(target_feat[t,:], [repetitions, 1]) - xs = ts + np.random.normal(size=ts.shape) @ L.T - loglik = utils.multiple_logpdfs_vec_input_single_cov(xs, template_feat, logdet, Us).squeeze() + xs = ts + rng.normal(size=ts.shape) @ L.T + loglik = utils.multiple_logpdfs_vec_input_single_cov( + xs,template_feat, logdet, Us).squeeze() logpost = loglik + np.log(prior) # normalise in log space for numerical stability logpost = logpost - logsumexp(logpost, axis=1, keepdims=True) # add numerical precision to avoid underflow (i.e. prob = 0) # it also function as a negligible lapse rate - logpost = np.logaddexp(logpost, np.log(np.finfo(loglik.dtype).eps)) - # normalise again (there is a better way but this is ok for now) + logpost = np.logaddexp( + logpost, np.log(np.finfo(loglik.dtype).eps)) + # normalise again + # (there is a better way but this is ok for now) logpost = logpost - logsumexp(logpost, axis=1, keepdims=True) posterior[t, :,:] = logpost else: for t in range(target_num): # for ta in range(target_num): # AWGN NOISE - x = np.random.multivariate_normal(target_feat[t,:], sigma) + x = rng.multivariate_normal(target_feat[t,:], sigma) # COMPUTE POSTERIOR # using vectorised solution - loglik = utils.multiple_logpdfs_vec_input_single_cov(np.expand_dims(x, axis=0), template_feat, logdet, Us).squeeze() + loglik = utils.multiple_logpdfs_vec_input_single_cov( + np.expand_dims(x, axis=0), + template_feat, + logdet, + Us, + ).squeeze() # post = np.exp(loglik+np.log(prior)) logpost = loglik + np.log(prior) # normalise logpost = logpost - logsumexp(logpost) # add numerical precision to avoid underflow (i.e. prob = 0) # it also function as a negligible lapse rate - logpost = np.logaddexp(logpost, np.log(np.finfo(loglik.dtype).eps)) - # normalise again (there is a better way but this is ok for now) + logpost = np.logaddexp(logpost, + np.log(np.finfo(loglik.dtype).eps)) + # normalise again + # (there is a better way but this is ok for now) logpost = logpost - logsumexp(logpost) # the solution above is faster than the for loop below but I am # keeping it for future reference and debugging # post = np.zeros(template_num) # for tp in range(template_num): - # # post[tp] = multivariate_normal.pdf(x, mean=template_feat[tp], cov=sigma) * prior[tp] + # # post[tp] = multivariate_normal.pdf( + # # x,mean=template_feat[tp], cov=sigma) * prior[tp] # # doing this speeds up stuff # u_diff = (x-template_feat[tp]) - # post[tp] = (np.exp(-0.5*u_diff @ sigma_inv @ u_diff.T))*prior[tp] + # post[tp] = ( + # np.exp(-0.5*u_diff @ sigma_inv @ u_diff.T))*prior[tp] # post /= np.sum(post) # logpost = np.log(post+np.finfo(post.dtype).eps) @@ -330,7 +407,8 @@ def estimate(self, posterior, sigma_motor=None): Returns ------- estimations : ndarray - Estimated directions in Cartesian coordinates (trials x repetitions x 3) + Estimated directions in Cartesian coordinates + (trials x repetitions x 3) """ repetitions = np.size(posterior, 1) trials = np.size(posterior, 0) @@ -350,7 +428,8 @@ def estimate(self, posterior, sigma_motor=None): if sigma_motor not in [False, 0]: for rt in range(repetitions): - estimations[:, rt, :] = utils.scatter_von_mises(estimations[:, rt, :], sigma_motor) + estimations[:, rt, :] = utils.scatter_von_mises( + estimations[:, rt, :], sigma_motor) return estimations @@ -390,7 +469,9 @@ def plot_cues(self, title='', fig=None, ax=None, clim=None, elev_min=None): else: im.set_clim(np.min(amps), np.max(amps)) - ax.set_xticks(np.interp([100, 1e3, 5e3, 1e4], self.freqs, np.arange(len(self.freqs)))) + ax.set_xticks(np.interp([100, 1e3, 5e3, 1e4], + self.freqs, + np.arange(len(self.freqs)))) ax.set_xticklabels([f'{freq:.0f}' for freq in [100, 1e3, 5e3, 1e4]]) ax.set_yticks(np.arange(len(elevations))) ax.set_yticklabels([f'{elev:.0f}' for elev in elevations]) @@ -399,7 +480,9 @@ def plot_cues(self, title='', fig=None, ax=None, clim=None, elev_min=None): def plot_post(self, posterior, estimations): amps = posterior.squeeze() - self.template.coords.plot(np.maximum(amps, np.log(np.finfo(amps.dtype).eps)), estimations.squeeze()) + self.template.coords.plot(np.maximum(amps, + np.log(np.finfo(amps.dtype).eps)), + estimations.squeeze()) def test_interp(): # doing it from scratch @@ -423,8 +506,14 @@ def test_interp(): fig, axes = plt.subplots(1, 2, figsize=(16, 6)) # Plot target and template features side by side with same color limits - am.plot_cues('- Target features (i.e. original)', fig=fig, ax=axes[0], clim=clim) - am.template.plot_cues('- Template features (i.e. interpolated)', fig=fig, ax=axes[1], clim=clim, elev_min=-45) + am.plot_cues('- Target features (i.e. original)', + fig=fig, + ax=axes[0], + clim=clim) + am.template.plot_cues('- Template features (i.e. interpolated)', + fig=fig, ax=axes[1], + clim=clim, + elev_min=-45) plt.tight_layout() plt.show() From 19fd618ee9b833e8edf2a27dc4a1a2e2a2fa2922 Mon Sep 17 00:00:00 2001 From: Emanuele Zanoni Date: Tue, 13 Jan 2026 19:45:37 +0100 Subject: [PATCH 2/6] fix(bayesian_listener): apply Ruff fixes (except D100, F841, FIX002 warnings) --- bayesian_listener/coordinates.py | 353 ++++++++++++++++++++++++++----- 1 file changed, 299 insertions(+), 54 deletions(-) diff --git a/bayesian_listener/coordinates.py b/bayesian_listener/coordinates.py index e5f4ed5..7867aec 100644 --- a/bayesian_listener/coordinates.py +++ b/bayesian_listener/coordinates.py @@ -17,12 +17,17 @@ class Coordinates: - lateral: [-90°, 90°] or [-π/2, π/2] rad - polar: [-90°, 270°) or [-π/2, 3π/2) rad - Note: The constructor automatically wraps angles to the correct ranges for each - coordinate system. SOFA files that use azimuth in [0°, 360°] will be wrapped - to [-180°, 180°]. + Note: The constructor automatically wraps angles to the correct ranges + for each coordinate system. + SOFA files that use azimuth in [0°, 360°] will be wrapped to [-180°, 180°]. """ - def __init__(self, sofa_file = None, positions = [], convention = 'cartesian', units = 'rad'): - if sofa_file != None: + + def __init__(self, + sofa_file = None, + positions = [], + convention = 'cartesian', + units = 'rad'): + if sofa_file is not None: # handle sofa input if isinstance(sofa_file, str): self.sofa_file = sofa_file @@ -41,8 +46,13 @@ def __init__(self, sofa_file = None, positions = [], convention = 'cartesian', u self.positions[:, 2] = 1.0 # force distance to be unitary self.convention = 'spherical' else: - if convention not in ['cartesian', 'spherical', 'horizontal-polar']: - raise ValueError(f'Specified "convention " is not supported: {convention}') + if convention not in [ + 'cartesian', + 'spherical', + 'horizontal-polar', + ]: + raise ValueError( + f'Specified "convention " is not supported: {convention}') # Convert to numpy array if needed positions = np.asarray(positions) @@ -54,15 +64,21 @@ def __init__(self, sofa_file = None, positions = [], convention = 'cartesian', u # Single position as 1D array, reshape to (1, 3) positions = positions.reshape(1, 3) else: - raise ValueError(f"1D positions array must have exactly 3 elements, got {positions.shape[0]}") + raise ValueError( + f"1D positions array must have exactly " + f"3 elements, got {positions.shape[0]}") elif positions.ndim == 2: if positions.shape[1] != 3: - raise ValueError(f"positions must be an Nx3 array, got shape {positions.shape}") + raise ValueError( + f"positions must be an Nx3 array, " + f"got shape {positions.shape}") else: - raise ValueError(f"positions must be 1D (size 3) or 2D (Nx3), got {positions.ndim}D array") + raise ValueError( + f"positions must be 1D (size 3) or 2D (Nx3), " + f"got {positions.ndim}D array") if (units == 'rad') & (np.abs(positions) > 2. * np.pi).any(): - warnings.warn('Coordinates ask for radiants!') + warnings.warn('Coordinates ask for radiants!', stacklevel=2) if units == 'deg': positions[:, :-1] = np.deg2rad(positions[:, :-1]) @@ -71,29 +87,77 @@ def __init__(self, sofa_file = None, positions = [], convention = 'cartesian', u self.convention = convention # Wrap angles to correct range for the coordinate system - if self.convention in ['spherical', 'horizontal-polar'] and np.size(self.positions) > 0: - self.positions = Coordinates._wrap_angles(self.positions, self.convention ) + if ( + self.convention in ['spherical', 'horizontal-polar'] + and np.size(self.positions) > 0 + ): + self.positions = Coordinates._wrap_angles(self.positions, + self.convention) + @classmethod def with_repetitions(cls, estimations, convention = 'cartesian'): + """ + Create a Coordinates instance from estimations with repetitions. + + Parameters + ---------- + estimations : np.ndarray + A 3D array of shape (N, reps, 3) containing N positions, + each with 'reps' repetitions of 3 coordinates. + convention : str + The coordinate system convention of the estimations. + Options: 'cartesian', 'spherical', 'horizontal-polar'. + + Returns + ------- + Coordinates : Coordinates + A Coordinates instance with the collapsed estimations. + """ # Collaps the repetitions dimension (from (N, reps, 3) to (N*reps, 3)) - assert estimations.ndim == 3, "Estimations should be a 3D array with shape (N, reps, 3)" - assert estimations.shape[2] == 3, "Estimations should have 3 coordinates in the last dimension" + assert estimations.ndim == 3, \ + "Estimations should be a 3D array with shape (N, reps, 3)" + assert estimations.shape[2] == 3, \ + "Estimations should have 3 coordinates in the last dimension" estimations = estimations.reshape(-1, 3) - # Create a new instance of Coordinates with the given positions and convention + # Create a new instance of Coordinates + # with the given positions and convention return cls(positions=estimations, convention=convention) def normalise(self): + """ + Normalize the positions by setting the radius to 1 + in spherical coordinates. + """ positions = self.convert('spherical') positions[:, 2] = 1 self.positions = positions self.convention = 'spherical' def convert(self, to_convention): - if to_convention not in ['cartesian', 'spherical', 'horizontal-polar']: - raise ValueError(f'Specified "to_convention" is not supported: {to_convention }') + """ + Convert the coordinates to a different convention. + + Parameters + ---------- + to_convention : str + The target coordinate system convention. + Options: 'cartesian', 'spherical', 'horizontal-polar'. + + Returns + ------- + np.ndarray : The converted coordinates. + """ + if to_convention not in [ + 'cartesian', + 'spherical', + 'horizontal-polar', + ]: + raise ValueError( + f"Specified 'to_convention' is not supported: " + f"{to_convention }") temp = self.positions.copy() if temp.ndim == 1: @@ -101,7 +165,9 @@ def convert(self, to_convention): if self.convention != to_convention: if self.convention in ['spherical']: - temp = Coordinates.sph2cart(temp[:, 0], temp[:, 1], temp[:, 2]) + temp = Coordinates.sph2cart(temp[:, 0], + temp[:, 1], + temp[:, 2]) elif self.convention == 'horizontal-polar': hor2sph_result = Coordinates.hor2sph(temp[:, 0], temp[:, 1]) az = hor2sph_result[:, 0] @@ -120,33 +186,91 @@ def convert(self, to_convention): return temp def sph(self): + """ + Get the spherical coordinates in degrees. + + Returns + ------- + np.ndarray : The spherical coordinates + (azimuth, elevation) in degrees. + """ positions = self.convert('spherical') positions[:, 0] = np.mod(positions[:, 0] + np.pi, 2 * np.pi) - np.pi return np.rad2deg(positions[:, (0, 1)]) def hpo(self): + """ + Get the horizontal-polar coordinates in degrees. + + Returns + ------- + np.ndarray : The horizontal-polar coordinates + (latitude, polar angle) in degrees. + """ positions = self.convert('horizontal-polar') positions[:, 0] = np.mod(positions[:, 0] + np.pi, 2 * np.pi) - np.pi return np.rad2deg(positions[:, (0, 1)]) def az(self): + """ + Get the azimuth angles in radians. + + Returns + ------- + np.ndarray : The azimuth angles in radians. + """ dirs = self.sph() return dirs[:, 0] def el(self): + """ + Get the elevation angles in radians. + + Returns + ------- + np.ndarray : The elevation angles in radians. + """ dirs = self.sph() return dirs[:, 1] def lat(self): + """ + Get the latitude angles in radians. + + Returns + ------- + np.ndarray : The latitude angles in radians. + """ dirs = self.hpo() return dirs[:, 0] def pol(self): + """ + Get the polar angles in radians. + + Returns + ------- + np.ndarray : The polar angles in radians. + """ dirs = self.hpo() return dirs[:, 1] - # Find the closest coordinates in the current set of coordinates def find(self, coords_search): + """ + Find the closest coordinates in the current set of coordinates. + + Parameters + ---------- + coords_search : Coordinates + The coordinates to search for. + + Returns + ------- + coords_found : Coordinates + The closest coordinates found. + idx : np.ndarray + The indices of the closest coordinates found. + """ pos = self.convert('cartesian') pos_search = coords_search.convert('cartesian') @@ -156,10 +280,23 @@ def find(self, coords_search): dist = np.sum((pos - pos_search[ii, :]) ** 2, axis=1) idx[ii] = np.argmin(dist) - coords_found = Coordinates(positions = self.positions[idx, :], convention = self.convention) + coords_found = Coordinates(positions = self.positions[idx, :], + convention = self.convention) return coords_found, idx def plot(self, values = [], points = None): + """ + Visualize the positions as a 3D scatter plot. + + Parameters + ---------- + values : array-like, optional + Numerical values used to color the scatter points. + If None or empty, all points are plotted with the same color. + points : array-like of shape (3,), optional + Optional 3D point to plot as an estimated direction + from the origin. + """ if values is None or len(values) == 0: values = np.ones(self.positions.shape[0]) else: @@ -168,11 +305,23 @@ def plot(self, values = [], points = None): r = self.convert('cartesian') fig = plt.figure() ax = fig.add_subplot(111, projection='3d') - ax.scatter(r[:,0], r[:, 1], r[:, 2], c = values, s=20, alpha = .5, label='Log posterior') + ax.scatter(r[:,0], + r[:, 1], + r[:, 2], + c = values, + s=20, + alpha = .5, + label='Log posterior', + ) ax.plot([0, 1], [0,0],zs=[0,0], c='red', label='Front direction') if points is not None: - ax.plot([0, points[0]], [0,points[1]],zs=[0,points[2]], c='blue', label='Estimated direction') + ax.plot(xs=[0, points[0]], + ys=[0, points[1]], + zs=[0, points[2]], + c='blue', + label='Estimated direction', + ) ax.view_init(elev=0, azim=0) ax.set_box_aspect([1, 1, 1]) @@ -182,6 +331,7 @@ def plot(self, values = [], points = None): plt.show() def print(self): + """Print the spherical coordinates in a grid format.""" grid = self.sph() print("Positions (spherical):") @@ -192,8 +342,33 @@ def print(self): def localization_error(self, estimations, metric, auxiliary_output=False): + """ + Compute the localization error between the current coordinates + and the provided estimations using the specified metric. + + Parameters + ---------- + estimations : Coordinates + The estimated coordinates to compare against. + metric : str or callable + The metric to use for error computation. + If a string, it should be a registered metric name. + If callable, it should be a function that takes two + Coordinates instances and returns the error value. + auxiliary_output : bool, optional + If True, return auxiliary output from the metric function. + Default is False. + + Returns + ------- + float or tuple : + The computed localization error. + If auxiliary_output is True, returns a tuple + (error_value, auxiliary_data). + """ if not isinstance(estimations, Coordinates): - raise ValueError("estimations must be an instance of Coordinates class") + raise ValueError( + "estimations must be an instance of Coordinates class") if self.positions.shape != estimations.positions.shape: raise ValueError("Shape mismatch") @@ -203,29 +378,41 @@ def localization_error(self, estimations, metric, auxiliary_output=False): # Case 2: metric is a string, but not registered in METRIC_FUNCTIONS if metric not in mt.METRIC_FUNCTIONS: - raise ValueError(f"Unknown metric: {metric}. Available metrics are: {list(mt.METRIC_FUNCTIONS.keys())}") + raise ValueError( + f"Unknown metric: {metric}. Available metrics are: " + f"{list(mt.METRIC_FUNCTIONS.keys())}") # Case 3: metric is a string and registered in METRIC_FUNCTIONS - expected_coord_convention = mt.get_metric_metadata(metric)['coord_convention'] + expected_coord_convention = \ + mt.get_metric_metadata(metric)['coord_convention'] expected_unit = mt.get_metric_metadata(metric)['input_unit'] - # Convert both self and estimations to the expected coordinate convention + # Convert both self and estimations + # to the expected coordinate convention converted_self = self.convert(expected_coord_convention) converted_estim = estimations.convert(expected_coord_convention) # TODO: Evaluate the unit if necessary - value, aux_out = mt.METRIC_FUNCTIONS[metric](converted_self, converted_estim) + value, aux_out = \ + mt.METRIC_FUNCTIONS[metric](converted_self, converted_estim) return (value, aux_out) if auxiliary_output else value @staticmethod def help_on_metric(name=None): + """ + Print help information about a specific metric + or list all available metrics if no name is provided. + """ mt.describe_metrics(name) @staticmethod def get_metric_metadata(name): + """ + Get metadata for a specific metric. + """ return mt.get_metric_metadata(name) @@ -234,32 +421,56 @@ def _wrap_angles(positions, coord_convention ): """ Wrap angles to the correct range for the specified coordinate system. - Parameters: - positions : np.ndarray - Array of shape (N, 3) containing positions in the specified coordinate system - coord_convention : str - convention of coordinate system: 'spherical' or 'horizontal-polar' - - Returns: - np.ndarray : positions with angles wrapped to the correct ranges + Parameters + ---------- + positions : np.ndarray + Array of shape (N, 3) containing positions in + the specified coordinate system + coord_convention : str + convention of coordinate system: + 'spherical' or 'horizontal-polar' + + Returns + ------- + np.ndarray : positions with angles wrapped to the correct ranges """ positions = positions.copy() if coord_convention == 'spherical': # Spherical: azimuth in [-pi, pi], elevation in [-pi/2, pi/2] positions[:, 0] = np.mod(positions[:, 0] + np.pi, 2*np.pi) - np.pi - # Elevation is naturally constrained by arcsin in cart2sph, but clip for safety + # Elevation is naturally constrained + # by arcsin in cart2sph, but clip for safety positions[:, 1] = np.clip(positions[:, 1], -np.pi/2, np.pi/2) elif coord_convention == 'horizontal-polar': - # Horizontal-polar: lateral in [-pi/2, pi/2], polar in [-pi/2, 3*pi/2) + # Horizontal-polar: + # lateral in [-pi/2, pi/2], polar in [-pi/2, 3*pi/2) positions[:, 0] = np.clip(positions[:, 0], -np.pi/2, np.pi/2) - positions[:, 1] = np.mod(positions[:, 1] + np.pi/2, 2*np.pi) - np.pi/2 + positions[:, 1] = \ + np.mod(positions[:, 1] + np.pi/2, 2*np.pi) - np.pi/2 return positions @staticmethod def sph2cart(azimuth, elevation, r): + """ + Transform spherical to cartesian coordinates. + + Parameters + ---------- + azimuth : np.ndarray + azimuth (in radians) + elevation : np.ndarray + elevation (in radians) + r : np.ndarray + radius (in meters) + + Returns + ------- + x, y, z : np.ndarray + cartesian coordinates (in meters) + """ x = r * np.cos(elevation) * np.cos(azimuth) y = r * np.cos(elevation) * np.sin(azimuth) z = r * np.sin(elevation) @@ -267,6 +478,27 @@ def sph2cart(azimuth, elevation, r): @staticmethod def cart2sph(x, y, z): + """ + Transform cartesian to spherical coordinates. + + Parameters + ---------- + x : np.ndarray + x coordinate (in meters) + y : np.ndarray + y coordinate (in meters) + z : np.ndarray + z coordinate (in meters) + + Returns + ------- + azimuth : np.ndarray + azimuth (in radians) + elevation : np.ndarray + elevation (in radians) + r : np.ndarray + radius (in meters) + """ r = np.sqrt(x**2 + y**2 + z**2) azimuth = np.arctan2(y, x) elevation = np.atan2(z, np.sqrt(x**2 + y**2)) @@ -277,13 +509,19 @@ def sph2hor(azi, ele): """ Transform spherical to horizontal-polar coordinates. - Parameters: - azi : azimuth (in radians) - ele : elevation (in radians) - - Returns: - lat : lateral angle (-pi/2 <= lat <= pi/2) - pol : polar angle (-pi/2 <= pol < 3*pi/2) + Parameters + ---------- + azi : np.ndarray + azimuth (in radians) + ele : np.ndarray + elevation (in radians) + + Returns + ------- + lat : np.ndarray + lateral angle (-pi/2 <= lat <= pi/2) + pol : np.ndarray + polar angle (-pi/2 <= pol < 3*pi/2) """ # Convert spherical to cartesian x, y, z = Coordinates.sph2cart(azi, ele, np.ones_like(azi)).T @@ -293,7 +531,8 @@ def sph2hor(azi, ele): y[np.abs(y) < np.finfo(float).eps] = 0 z[np.abs(z) < np.finfo(float).eps] = 0 - # Interpret horizontal polar format as rotated spherical coordinates with negative azimuth direction + # Interpret horizontal polar format as rotated spherical coordinates + # with negative azimuth direction pol, nlat, r = Coordinates.cart2sph(x, z, -y).T lat = -nlat @@ -307,13 +546,19 @@ def hor2sph(lat, pol): """ Transform horizontal-polar to spherical coordinates. - Parameters: - lat : lateral angle (-pi/2 <= lat <= pi/2, in radians) - pol : polar angle (-pi/2 <= pol < 3*pi/2, in radians) - - Returns: - azi : azimuth (-pi <= azi <= pi, in radians) - ele : elevation (-pi/2 <= ele <= pi/2, in radians) + Parameters + ---------- + lat : np.ndarray + lateral angle (-pi/2 <= lat <= pi/2, in radians) + pol : np.ndarray + polar angle (-pi/2 <= pol < 3*pi/2, in radians) + + Returns + ------- + azi : np.ndarray + azimuth (-pi <= azi <= pi, in radians) + ele : np.ndarray + elevation (-pi/2 <= ele <= pi/2, in radians) """ x, nz, y = Coordinates.sph2cart(-pol, lat, np.ones_like(lat)).T From 21ccfdd130d18466c55ce1c70b369c4561df85e4 Mon Sep 17 00:00:00 2001 From: Emanuele Zanoni Date: Tue, 13 Jan 2026 19:59:45 +0100 Subject: [PATCH 3/6] fix(metrics): apply Ruff fixes (except D100 warning) --- bayesian_listener/metrics.py | 123 +++++++++++++++++++++++++++++------ 1 file changed, 102 insertions(+), 21 deletions(-) diff --git a/bayesian_listener/metrics.py b/bayesian_listener/metrics.py index 4311c7d..273ca8f 100644 --- a/bayesian_listener/metrics.py +++ b/bayesian_listener/metrics.py @@ -3,23 +3,59 @@ # Shared dictionary to hold metric functions and their metadata METRIC_FUNCTIONS = {} -def register_metric(name, coord_convention, input_unit, output_unit=None, description=None, **extra_metadata): +def register_metric(name, + coord_convention, + input_unit, + output_unit=None, + description=None, + **extra_metadata, + ): + """ + Decorator to register a metric function with metadata. + + Parameters + ---------- + name : str + Name of the metric. + coord_convention : str + Coordinate convention used (e.g., 'horizontal-polar'). + input_unit : str + Unit of the input data (e.g., 'radians'). + output_unit : str, optional + Unit of the output data (e.g., 'radians', 'percentage'). + description : str, optional + Description of the metric. + **extra_metadata : dict + Additional metadata to store. + + Returns + ------- + decorator : function + Decorator that registers the metric function. + """ def decorator(func): + """ + Decorator that registers the metric function with metadata. + """ def wrapped(*args, **kwargs): + """ + Wrapper to ensure uniform output format. + """ result = func(*args, **kwargs) if isinstance(result, tuple): value, auxiliary_output = result else: value = result auxiliary_output = {} - return value, auxiliary_output # Every function is uniformly formatted to return a tuple + # Every function is uniformly formatted to return a tuple + return value, auxiliary_output wrapped._metadata = { 'name': name, 'coord_convention': coord_convention, 'input_unit': input_unit, 'output_unit': output_unit, 'description': description, - **extra_metadata + **extra_metadata, } METRIC_FUNCTIONS[name] = wrapped return wrapped @@ -27,13 +63,35 @@ def wrapped(*args, **kwargs): def get_metric_metadata(name): + """ + Retrieve metadata for a registered metric. + + Parameters + ---------- + name : str + Name of the metric. + + Returns + ------- + metadata : dict + Metadata dictionary for the metric. + """ func = METRIC_FUNCTIONS.get(name) if func is None: raise ValueError(f"Metric '{name}' not found.") - return func._metadata.copy() # Return a copy to prevent external modification + # Return a copy to prevent external modification + return func._metadata.copy() def describe_metrics(name=None): + """ + Print descriptions of registered metrics. + + Parameters + ---------- + name : str, optional + Name of the metric to describe. If None, lists all metrics. + """ if name: info = get_metric_metadata(name) print(f"Metric: {name}") @@ -43,7 +101,8 @@ def describe_metrics(name=None): print("Available metrics:") for name in METRIC_FUNCTIONS.keys(): print(f" {name}: {get_metric_metadata(name)['description']}") - print("Use describe_metrics(name) to get details for a specific metric.") + print( + "Use describe_metrics(name) to get details for a specific metric.") def wrap_to_pi(rad): @@ -53,14 +112,14 @@ def wrap_to_pi(rad): def wrap_polar_angle(angle_rad): """ - Wrap polar angles to the range [-π/2, 3π/2) ≡ [-90°, 270°) + Wrap polar angles to the range [-π/2, 3π/2) ≡ [-90°, 270°). """ return (angle_rad + np.pi / 2) % (2 * np.pi) - np.pi / 2 -# ------------------------------------------------------------------------ +# ----------------------------------------------------------------------------- # Metric Functions @@ -71,12 +130,18 @@ def wrap_polar_angle(angle_rad): output_unit="radians", description=( "Lateral RMS error (in radians).\n\t" - "RMS of the difference between response and target lateral angles within ±60° lateral.\n\t" + "RMS of the difference between response and target lateral angles\n\t" + "within ±60° lateral.\n\t" "See rms lateral error in Middlebrooks (1999)"), - ylabel="Lateral RMS error (rad)" + ylabel="Lateral RMS error (rad)", ) def rmsL(true, est): - lat_true = wrap_to_pi(true[..., 0]) # lateral in [-π, π), then restrict to [-π/2, π/2] + """ + Compute lateral RMS error within ±60° lateral. + More details in the decorator above. + """ + # lateral in [-π, π), then restrict to [-π/2, π/2] + lat_true = wrap_to_pi(true[..., 0]) lat_true = np.clip(lat_true, -np.pi/2, np.pi/2) # enforce [-π/2, π/2] lat_est = wrap_to_pi(est[..., 0]) @@ -97,23 +162,31 @@ def rmsL(true, est): output_unit="radians", description=( "RMS polar error (local, central responses only).\n\t" - "Root mean square of polar angle error, restricted to responses with:\n\t" + "Root mean square of polar angle error,\n\t" + "restricted to responses with:\n\t" "- lateral response within ±30° (±π/6 radians)\n\t" "- polar error less than 90° (π/2 radians).\n\t" "Based on definition in Middlebrooks (1999)." ), - ylabel="Local central RMS polar error (rad)" + ylabel="Local central RMS polar error (rad)", ) def rmsPmedianlocal(true, est): - lat_est = wrap_to_pi(est[..., 0]) # lateral in [-π, π), then restrict to [-π/2, π/2] - assert np.all(np.abs(lat_est) <= np.pi/2), "Lateral angles must be in [-π/2, π/2]" + """ + Compute local RMS polar error within ±30° lateral and polar error < 90°. + More details in the decorator above. + """ + # lateral in [-π, π), then restrict to [-π/2, π/2] + lat_est = wrap_to_pi(est[..., 0]) + assert np.all(np.abs(lat_est) <= np.pi/2), \ + "Lateral angles must be in [-π/2, π/2]" pol_true = wrap_polar_angle(true[..., 1]) # polar in [-π/2, 3π/2) pol_est = wrap_polar_angle(est[..., 1]) # 1. Select central responses: lateral response within ±30° central_mask = np.abs(lat_est) <= np.deg2rad(30) - assert np.any(central_mask), "No central responses found within ±30° lateral range." + assert np.any(central_mask), \ + "No central responses found within ±30° lateral range." # 2. Exclude responses with polar error greater than 90° polar_diff = wrap_to_pi(pol_est - pol_true)[central_mask] @@ -138,19 +211,27 @@ def rmsPmedianlocal(true, est): ylabel="Quadrant errors (%)", auxiliary_output={ 'confusion_count': 'Number of confusions (polar error ≥ 90°)', - 'response_count': 'Number of responses within the lateral range (|lat| ≤ 30°)' - } + 'response_count': \ + 'Number of responses within the lateral range (|lat| ≤ 30°)', + }, ) def querrMiddlebrooks(true, est): - lat_est = wrap_to_pi(est[..., 0]) # lateral in [-π, π), then restrict to [-π/2, π/2] - assert np.all(np.abs(lat_est) <= np.pi/2), "Lateral angles must be in [-π/2, π/2]" + """ + Compute quadrant error rate as defined in Middlebrooks (1999). + More details in the decorator above. + """ + # lateral in [-π, π), then restrict to [-π/2, π/2] + lat_est = wrap_to_pi(est[..., 0]) + assert np.all(np.abs(lat_est) <= np.pi/2), \ + "Lateral angles must be in [-π/2, π/2]" pol_true = wrap_polar_angle(true[..., 1]) # polar in [-π/2, 3π/2) pol_est = wrap_polar_angle(est[..., 1]) # 1. Filter central responses: lateral response within ±30° central_mask = np.abs(lat_est) <= np.deg2rad(30) - assert np.any(central_mask), "No central responses found within ±30° lateral range." + assert np.any(central_mask), \ + "No central responses found within ±30° lateral range." # 2. Compute polar error and count confusions (polar error ≥ 90°) polar_error = np.abs(wrap_to_pi(pol_est - pol_true))[central_mask] @@ -158,4 +239,4 @@ def querrMiddlebrooks(true, est): n_total = np.int64(len(polar_error)) qerr = 100 * n_confusions / n_total - return qerr, {'confusion_count': n_confusions, 'response_count': n_total} \ No newline at end of file + return qerr, {'confusion_count': n_confusions, 'response_count': n_total} From ce5af7e0d05f81d0cdc5a5451783b24007b046ba Mon Sep 17 00:00:00 2001 From: Emanuele Zanoni Date: Tue, 13 Jan 2026 20:13:20 +0100 Subject: [PATCH 4/6] fix(resample): apply Ruff fixes (except ERA and D warnings) --- bayesian_listener/resample.py | 124 ++++++++++++++++++++++++---------- 1 file changed, 88 insertions(+), 36 deletions(-) diff --git a/bayesian_listener/resample.py b/bayesian_listener/resample.py index 755763f..905be8e 100644 --- a/bayesian_listener/resample.py +++ b/bayesian_listener/resample.py @@ -10,7 +10,10 @@ # ----------------------------------------------------------------------------- # as in barumerli2023 -def resample_barumerli2023(values, coords_in, dirs = None, flag_regularisation = True): +def resample_barumerli2023(values, + coords_in, + dirs = None, + flag_regularisation = True): """ Resample using spherical harmonics as in Barumerli et al. 2023. @@ -18,7 +21,8 @@ def resample_barumerli2023(values, coords_in, dirs = None, flag_regularisation = ---------- values : array or list of arrays Single cue array of shape (n_dirs, ...) or list of cue arrays. - If list, each array must have shape (n_dirs, ...) where first dimension matches. + If list, each array must have shape (n_dirs, ...) + where first dimension matches. coords_in : Coordinates Source coordinates dirs : array or Coordinates, optional @@ -45,7 +49,8 @@ def resample_barumerli2023(values, coords_in, dirs = None, flag_regularisation = if dirs is None: # %% Generate t-design points - # the advantage of using this is that the weights are equal when integrating + # the advantage of using this is that + # the weights are equal when integrating dirs = spaudiopy.grids.load_n_design(64)# 2112 equally distant points # assert(N_SH < N_dirs, ... @@ -68,7 +73,9 @@ def resample_barumerli2023(values, coords_in, dirs = None, flag_regularisation = zen = np.pi - ele # get SH basis on new directions - int_new = spaudiopy.sph.sh_matrix(N_sph, dirs_SH[:, 0], dirs_SH[:, 1], sh_type='real') + int_new = spaudiopy.sph.sh_matrix(N_sph, dirs_SH[:, 0], + dirs_SH[:, 1], + sh_type='real') # Ensure all cues are at least 2D cues = [c[:, np.newaxis] if c.ndim == 1 else c for c in cues] @@ -85,7 +92,9 @@ def resample_barumerli2023(values, coords_in, dirs = None, flag_regularisation = # get SH basis on old directions Y_N_tik = spaudiopy.sph.sh_matrix(N_sph, azi, zen, 'real') # Compute regularized inverse once - Y_inv_reg = np.linalg.solve(np.transpose(Y_N_tik)@Y_N_tik+lambda_val*SIG, np.transpose(Y_N_tik)) + Y_inv_reg = np.linalg.solve( + np.transpose(Y_N_tik)@Y_N_tik+lambda_val*SIG, + np.transpose(Y_N_tik)) # Transform all cues to SH domain cues_SH = [Y_inv_reg @ c for c in cues] @@ -116,17 +125,21 @@ def build_Y(dirs, N): return Y def build_bau_damping(N): - """Bau et al. damping: D_ii = 1 + n(n+1)""" + """Bau et al. damping: D_ii = 1 + n(n+1).""" num_coeffs = (N + 1) ** 2 D = np.zeros((num_coeffs, num_coeffs)) idx = 0 for n in range(N + 1): - for m in range(-n, n + 1): + for _ in range(-n, n + 1): D[idx, idx] = 1 + n * (n + 1) idx += 1 return D -def find_max_order(dirs, thresh=12.25, N_max=35, regularised=True, epsilon = 1e-2): +def find_max_order(dirs, + thresh=12.25, + N_max=35, + regularised=True, + epsilon=1e-2): """ Return the largest N ≤ N_max such that cond(Y(dirs,N)) < thresh. @@ -169,8 +182,8 @@ def solve_sh(Y, H): def interpolate_HRTF(query_dirs, C, N): """ - query_dirs: (Q,2) array of (az,el) - returns: (Q,F) matrix of interpolated HRTF magnitudes + query_dirs: (Q,2) array of (az,el). + returns: (Q,F) matrix of interpolated HRTF magnitudes. """ Yq = build_Y(query_dirs, N) return Yq @ C @@ -178,18 +191,19 @@ def interpolate_HRTF(query_dirs, C, N): # this is the first attempt to interpolate with SH def resample_old(values, coords_in, dirs = None, plot_grid=False): """ - Resample function with regularisation and control on condition number + Resample function with regularisation and control on condition number. :param values: Description :param coords_in: Description :param dirs: Description :param plot_grid: Description """ - N_sph = 15 + # N_sph = 15 if dirs is None: # %% Generate t-design points - # the advantage of using this is that the weights are equal when integrating + # the advantage of using this is that + # the weights are equal when integrating dirs = spaudiopy.grids.load_n_design(64)# 2112 equally distant points else: assert(isinstance(dirs, Coordinates)) @@ -203,7 +217,7 @@ def resample_old(values, coords_in, dirs = None, plot_grid=False): dirs_original = c[:, (0, 1)] # move to navigation coordinates - azi = dirs_original[:, 0] # must be in [0, 2*pi] + # azi = dirs_original[:, 0] # must be in [0, 2*pi] ele = dirs_original[:, 1] if np.any(ele < 0): ele = ele + np.pi/2 # must be in [0, pi]. @@ -302,7 +316,8 @@ def complement_sampling(coordinates): if min_elevation > 0: warnings.warn( 'Detected positive minimum elevation during resampling.' - 'Manual resampling might be required') + 'Manual resampling might be required', + stacklevel=2) # find and add mirror points with added 0.1 degree safety margin mask = coordinates.elevation > -min_elevation + 0.0017 @@ -344,7 +359,7 @@ def resample_two_step(cues, coordinates, template, second_step): # check input format if not isinstance(cues, (list, tuple)): - cues = [cues, ] + cues = [cues] passed_list = False else: passed_list = True @@ -352,11 +367,15 @@ def resample_two_step(cues, coordinates, template, second_step): # Convert custom Coordinates to pyfar.Coordinates if needed if isinstance(coordinates, Coordinates): coords_cart = coordinates.convert('cartesian') - coordinates = pf.Coordinates.from_cartesian(coords_cart[:, 0], coords_cart[:, 1], coords_cart[:, 2]) + coordinates = pf.Coordinates.from_cartesian(coords_cart[:, 0], + coords_cart[:, 1], + coords_cart[:, 2]) # Convert template to pyfar.Coordinates if it's a numpy array if isinstance(template, np.ndarray): - template = pf.Coordinates.from_cartesian(template[:, 0], template[:, 1], template[:, 2]) + template = pf.Coordinates.from_cartesian(template[:, 0], + template[:, 1], + template[:, 2]) # HRTF measurement grids usually lack the bottom, so we add it coordinates_complemented, mask = complement_sampling(coordinates) @@ -444,7 +463,8 @@ def resample(cues, coordinates, method='SH'): ---------- cues : array or list of arrays Single cue array of shape (n_dirs, ...) or list of cue arrays. - If list, each array must have shape (n_dirs, ...) where first dimension matches. + If list, each array must have shape (n_dirs, ...) + where first dimension matches. coordinates : Coordinates Source coordinates method : str @@ -454,7 +474,8 @@ def resample(cues, coordinates, method='SH'): ------- result : array or list of arrays Resampled cues. Returns same type (single array or list) as input. - Each array has shape (n_template_dirs, ...) matching input except first dimension. + Each array has shape (n_template_dirs, ...) + matching input except first dimension. template_coords : Coordinates Output coordinates """ @@ -462,23 +483,32 @@ def resample(cues, coordinates, method='SH'): if method.lower() == 'barycentric': result = resample_two_step(cues, coordinates, template, 'barycentric') - template_coords = Coordinates(positions=template, convention='cartesian') + template_coords = Coordinates(positions=template, + convention='cartesian') elif method.lower() == 'sh': result = resample_two_step(cues, coordinates, template, 'sh') - template_coords = Coordinates(positions=template, convention='cartesian') + template_coords = Coordinates(positions=template, + convention='cartesian') elif method.lower() == 'shmax': result = resample_two_step(cues, coordinates, template, 'shmax') - template_coords = Coordinates(positions=template, convention='cartesian') + template_coords = Coordinates(positions=template, + convention='cartesian') elif method.lower() == 'barumerli2023': - result, template_coords = resample_barumerli2023(cues, coordinates, template) + result, template_coords = resample_barumerli2023(cues, + coordinates, + template) else: raise ValueError(f"Unknown resample method: {method}") return result, template_coords -def plot_resampling_grid(coords_meas_cart, dirs_virt, missing_mask, z_min_meas): +def plot_resampling_grid(coords_meas_cart, + dirs_virt, + missing_mask, + z_min_meas): """ - Plot the resampling grid showing measured directions (black) and added directions (red). + Plot the resampling grid showing measured directions (black) + and added directions (red). Parameters ---------- @@ -499,13 +529,25 @@ def plot_resampling_grid(coords_meas_cart, dirs_virt, missing_mask, z_min_meas): ax1 = fig.add_subplot(121, projection='3d') # Plot measured directions (black) - ax1.scatter(coords_meas_cart[:, 0], coords_meas_cart[:, 1], coords_meas_cart[:, 2], - c='black', s=20, alpha=0.6, label=f'Measured (n={len(coords_meas_cart)})') + ax1.scatter(coords_meas_cart[:, 0], + coords_meas_cart[:, 1], + coords_meas_cart[:, 2], + c='black', + s=20, + alpha=0.6, + label=f'Measured (n={len(coords_meas_cart)})', + ) # Plot added directions (red) dirs_added = dirs_virt[missing_mask] - ax1.scatter(dirs_added[:, 0], dirs_added[:, 1], dirs_added[:, 2], - c='red', s=20, alpha=0.6, label=f'Added (n={np.sum(missing_mask)})') + ax1.scatter(dirs_added[:, 0], + dirs_added[:, 1], + dirs_added[:, 2], + c='red', + s=20, + alpha=0.6, + label=f'Added (n={np.sum(missing_mask)})', + ) # Draw horizontal plane at z_min_meas xx, yy = np.meshgrid(np.linspace(-1, 1, 10), np.linspace(-1, 1, 10)) @@ -524,10 +566,19 @@ def plot_resampling_grid(coords_meas_cart, dirs_virt, missing_mask, z_min_meas): cm = Coordinates(positions = coords_meas_cart, convention = 'cartesian') cd = Coordinates(positions = dirs_added, convention = 'cartesian') - ax2.scatter(cm.az(), cm.el(), - c='black', s=20, alpha=0.6, label=f'Measured (n={len(coords_meas_cart)})') - ax2.scatter(cd.az(), cd.el(), - c='red', s=20, alpha=0.6, label=f'Added (n={np.sum(missing_mask)})') + ax2.scatter(cm.az(), + cm.el(), + c='black', + s=20, + alpha=0.6, + label=f'Measured (n={len(coords_meas_cart)})', + ) + ax2.scatter(cd.az(), + cd.el(), + c='red', + s=20, alpha=0.6, + label=f'Added (n={np.sum(missing_mask)})', + ) ax2.set_xlabel('Lateral (deg)') ax2.set_ylabel('Polar (deg)') ax2.set_title('Top View: Resampling Grid') @@ -538,7 +589,8 @@ def plot_resampling_grid(coords_meas_cart, dirs_virt, missing_mask, z_min_meas): plt.tight_layout() plt.show() - print(f"Grid statistics:") + print("Grid statistics:") print(f" Measured directions: {len(coords_meas_cart)}") print(f" Added directions (z < {z_min_meas:.3f}): {np.sum(missing_mask)}") - print(f" Total directions for interpolation: {len(coords_meas_cart) + np.sum(missing_mask)}") + print(f" Total directions for interpolation: " + f"{len(coords_meas_cart) + np.sum(missing_mask)}") From ed818932058da01d55228f6123b60c41a8541681 Mon Sep 17 00:00:00 2001 From: Emanuele Zanoni Date: Tue, 13 Jan 2026 20:28:39 +0100 Subject: [PATCH 5/6] fix(utils): apply Ruff fixes (except ERA, D, F811 warnings) --- bayesian_listener/utils.py | 115 ++++++++++++++++++++++++------------- 1 file changed, 75 insertions(+), 40 deletions(-) diff --git a/bayesian_listener/utils.py b/bayesian_listener/utils.py index 14e3383..8bcbf2f 100644 --- a/bayesian_listener/utils.py +++ b/bayesian_listener/utils.py @@ -7,6 +7,8 @@ from pathlib import Path import pickle import datetime +import psutil +import os # feature functions def mag2db(mag): @@ -16,7 +18,8 @@ def erb_space(freq_range=[7e2, 18e3], erb_spacing=1): #translated from amt/audspacebw.m focusing on ERB-rate scale # Convert frequency limits to auditory scale (ERB-rate scale) - audlimits = 9.2645 * np.sign(freq_range) * np.log(1 + np.abs(freq_range) * 0.00437) + audlimits = \ + 9.2645 * np.sign(freq_range) * np.log(1 + np.abs(freq_range) * 0.00437) audrange = audlimits[1] - audlimits[0] # Calculate number of points (excluding final point) @@ -32,7 +35,8 @@ def erb_space(freq_range=[7e2, 18e3], erb_spacing=1): n += 1 # Convert auditory scale points back to Hz - fc = (1 / 0.00437) * np.sign(audpoints) * (np.exp(np.abs(audpoints) / 9.2645) - 1) + fc = (1 / 0.00437) * np.sign(audpoints) * \ + (np.exp(np.abs(audpoints) / 9.2645) - 1) return fc @@ -46,7 +50,7 @@ def gammatone( ): """ Complex-valued, all-pole gammatone filter coefficients as in Lyon, 1997. - The function has been taken from the AMT/gammatone.m + The function has been taken from the AMT/gammatone.m. Parameters ---------- @@ -62,8 +66,8 @@ def gammatone( scale : {"0dBforall", "6dBperoctave"} Amplitude scaling mode. phase : {"causalphase", "peakphase", "exppeakphase"} - Phase option. "exppeakphase" additionally aligns maxima using an impulse - response simulation (requires SciPy). + Phase option. "exppeakphase" additionally aligns maxima using an + impulse response simulation (requires SciPy). Returns ------- @@ -95,7 +99,8 @@ def gammatone( raise ValueError("scale must be '0dBforall' or '6dBperoctave'.") if phase not in {"causalphase", "peakphase", "exppeakphase"}: - raise ValueError("phase must be 'causalphase', 'peakphase', or 'exppeakphase'.") + raise ValueError( + "phase must be 'causalphase', 'peakphase', or 'exppeakphase'.") # ---- bandwidth multiplier (match the MATLAB formula literally) ---- if betamul is None: @@ -125,7 +130,7 @@ def gammatone( # ---- design each channel ---- for i in range(nch): # Complex pole location (all-pole, repeated n times) - atilde = np.exp(-2 * np.pi * beta[i] / fs - 1j * 2 * np.pi * fc[i] / fs) + atilde = np.exp(-2*np.pi * beta[i] / fs - 1j * 2*np.pi * fc[i] / fs) # Denominator from repeated root (length n+1, leading 1) # np.poly takes roots and returns monic polynomial coefficients. @@ -156,7 +161,8 @@ def gammatone( envmax = np.argmax(np.abs(tmp)) sigmax = np.argmax(tmp) # Equation analogous to the MATLAB code: - phi_delay = fc[i] * (-2 * np.pi - np.pi / 4.0) * (envmax - sigmax) / fs + phi_delay = \ + fc[i] * (-2*np.pi - np.pi / 4.0) * (envmax - sigmax) / fs b_i = b_i * np.exp(1j * phi_delay) # Store results @@ -172,11 +178,13 @@ def itdestimator(signals, fs=None): """ Estimate ITD from the given stimulus. - Parameters: + Parameters + ---------- Obj : 3D numpy array or object with IR data fs : Sampling rate (required if Obj is a 3D array) - Returns: + Returns + ------- toa_diff : Time of arrival difference """ @@ -221,7 +229,8 @@ def itdestimator(signals, fs=None): def scatter_von_mises(dirs, sigma_m): assert dirs.shape[1] == 3 or dirs.size == 3 - assert sigma_m >= 5, 'sensorimotor concentration too small and can lead to complex values' + assert sigma_m >= 5, \ + 'sensorimotor concentration too small and can lead to complex values' dirs = np.squeeze(dirs) @@ -241,7 +250,8 @@ def scatter_von_mises(dirs, sigma_m): return dirs_new def randvmf(kappa, mu, seed = None): - np.random.seed(seed) + rng = np.random.default_rng(seed) + assert mu is not None assert kappa > 0 @@ -255,11 +265,11 @@ def randvmf(kappa, mu, seed = None): # Rubinstein 81, p.39, Fisher 87, p.59 kappaS = np.sign(kappa) kappa = abs(kappa) - U = np.random.rand() + U = rng.rand() x = np.log(2. * U * np.sinh(kappa) + np.exp(-kappa)) / kappa x = kappaS * x - psi = 2. * np.pi * np.random.rand() + psi = 2. * np.pi * rng.rand() s_x = np.sqrt(1. - x**2.) y = np.array([np.cos(psi) * s_x, np.sin(psi) * s_x, x]) @@ -314,7 +324,8 @@ def multiple_logpdfs_vec_input(xs, means, covs): Us = vecs * np.sqrt(valsinvs)[:, None] devs = xs[:, None, :] - means[None, :, :] - # Use `einsum` for matrix-vector multiplications across the first dimension. + # Use `einsum` for matrix-vector multiplications + # across the first dimension. devUs = np.einsum('jnk,nki->jni', devs, Us) # Compute the Mahalanobis distance by squaring each term and summing. @@ -334,13 +345,16 @@ def multiple_logpdfs_vec_input_single_cov(xs, means, logdet, Us): https://gregorygundersen.com/blog/2020/12/12/group-multivariate-normal-pdf/ The big idea is to do one intensive operation, eigenvalue decomposition, - and then use that decomposition to compute the matrix inverse and determinant cheaply. + and then use that decomposition to compute the matrix inverse + and determinant cheaply. """ devs = xs[:, None, :] - means[None, :, :] - # Use `einsum` for matrix-vector multiplications across the first dimension. - # devUs = np.einsum('jnk,ki->jni', devs, Us) -> using this notation is very slow (twice as much) + # Use `einsum` for matrix-vector multiplications + # across the first dimension. + # devUs = np.einsum('jnk,ki->jni', devs, Us) -> + # using this notation is very slow (twice as much) devUs = devs @ Us # Compute the Mahalanobis distance by squaring each term and summing. @@ -361,7 +375,8 @@ def multiple_logpdfs_vec_input_single_cov(xs, means, logdet, Us): https://gregorygundersen.com/blog/2020/12/12/group-multivariate-normal-pdf/ The big idea is to do one intensive operation, eigenvalue decomposition, - and then use that decomposition to compute the matrix inverse and determinant cheaply. + and then use that decomposition to compute the matrix inverse + and determinant cheaply. """ n_samples = xs.shape[0] @@ -394,13 +409,17 @@ def multiple_logpdfs_vec_input_single_cov(xs, means, logdet, Us): return out @jit(nopython=True, parallel=True) -def multiple_logpdfs_vec_input_single_cov_diagonal(xs, means, logdet, sigma_inv_diag): +def multiple_logpdfs_vec_input_single_cov_diagonal(xs, + means, + logdet, + sigma_inv_diag): """ Optimized version for diagonal covariance. sigma_inv_diag: (P,) array of 1/sigma_i² """ devs = xs[:, None, :] - means[None, :, :] # (n_samples, n_means, P) - # Mahalanobis distance for diagonal cov: dev.T @ inv(Σ) @ dev = sum(dev² / sigma²) + # Mahalanobis distance for diagonal cov: + # dev.T @ inv(Σ) @ dev = sum(dev² / sigma²) mahas = np.sum(devs**2 * sigma_inv_diag, axis=2) # (n_samples, n_means) dim = xs.shape[1] @@ -468,7 +487,10 @@ def clear_cache(sofa_file=None): else: print("✓ No cache index found") -def load_from_cache(cache_dir, sofa_file, attributes_to_restore, interpolation='SH'): +def load_from_cache(cache_dir, + sofa_file, + attributes_to_restore, + interpolation='SH'): """ Try to load cached data from pickle file. @@ -502,7 +524,8 @@ def load_from_cache(cache_dir, sofa_file, attributes_to_restore, interpolation=' # Handle backward compatibility: add 'interpolation' column if missing if 'interpolation' not in cache_df.columns: - cache_df['interpolation'] = 'SH' # Assume old caches used default SH method + # Assume old caches used default SH method + cache_df['interpolation'] = 'SH' # Check for matching entry (including interpolation method) match = cache_df[(cache_df['sofa_name'] == sofa_name) & @@ -522,12 +545,13 @@ def load_from_cache(cache_dir, sofa_file, attributes_to_restore, interpolation=' cached_data = pickle.load(f) # Validate that all required attributes are present - missing = [attr for attr in attributes_to_restore if attr not in cached_data] + missing = [ + attr for attr in attributes_to_restore if attr not in cached_data] if missing: print(f"⚠ Cache missing attributes: {missing}") return None - print(f"✓ Cache loaded successfully") + print("✓ Cache loaded successfully") return cached_data except Exception as e: @@ -567,9 +591,16 @@ def save_to_cache(cache_dir, sofa_file, data_to_cache, interpolation='SH'): cache_df = pd.read_csv(cache_index_file) # Handle backward compatibility: add 'interpolation' column if missing if 'interpolation' not in cache_df.columns: - cache_df['interpolation'] = 'SH' # Assume old caches used default SH method + # Assume old caches used default SH method + cache_df['interpolation'] = 'SH' else: - cache_df = pd.DataFrame(columns=['sofa_name', 'file_hash', 'interpolation', 'pkl_file', 'timestamp']) + cache_df = pd.DataFrame(columns=[ + 'sofa_name', + 'file_hash', + 'interpolation', + 'pkl_file', + 'timestamp', + ]) # Check if entry already exists (same name, hash, and interpolation) match = cache_df[(cache_df['sofa_name'] == sofa_name) & @@ -577,7 +608,11 @@ def save_to_cache(cache_dir, sofa_file, data_to_cache, interpolation='SH'): (cache_df['interpolation'] == interpolation)] # Create new pickle filename with interpolation method and timestamp - pkl_filename = f"{Path(sofa_name).stem}_{file_hash[:8]}_{interpolation}_{timestamp}.pkl" + pkl_filename = ( + f"{Path(sofa_name).stem}_{file_hash[:8]}" + f"_{interpolation}_{timestamp}.pkl" + ) + pkl_path = cache_dir / pkl_filename if not match.empty: @@ -601,10 +636,12 @@ def save_to_cache(cache_dir, sofa_file, data_to_cache, interpolation='SH'): cache_df.loc[match.index[0], 'pkl_file'] = pkl_filename cache_df.loc[match.index[0], 'timestamp'] = timestamp else: - # Remove old entries for same file with different hash (but keep different interpolations) - old_entries = cache_df[(cache_df['sofa_name'] == sofa_name) & - (cache_df['file_hash'] != file_hash) & - (cache_df['interpolation'] == interpolation)] + # Remove old entries for same file with different hash + # (but keep different interpolations) + old_entries = cache_df[ + (cache_df['sofa_name'] == sofa_name) & + (cache_df['file_hash'] != file_hash) & + (cache_df['interpolation'] == interpolation)] # Delete old pickle files with different hashes for _, row in old_entries.iterrows(): @@ -614,9 +651,10 @@ def save_to_cache(cache_dir, sofa_file, data_to_cache, interpolation='SH'): print(f"→ Removed outdated cache: {row['pkl_file']}") # Remove old entries from dataframe - cache_df = cache_df[~((cache_df['sofa_name'] == sofa_name) & - (cache_df['file_hash'] != file_hash) & - (cache_df['interpolation'] == interpolation))] + cache_df = cache_df[~( + (cache_df['sofa_name'] == sofa_name) & + (cache_df['file_hash'] != file_hash) & + (cache_df['interpolation'] == interpolation))] # Add new entry new_row = pd.DataFrame([{ @@ -624,12 +662,12 @@ def save_to_cache(cache_dir, sofa_file, data_to_cache, interpolation='SH'): 'file_hash': file_hash, 'interpolation': interpolation, 'pkl_file': pkl_filename, - 'timestamp': timestamp + 'timestamp': timestamp, }]) cache_df = pd.concat([cache_df, new_row], ignore_index=True) cache_df.to_csv(cache_index_file, index=False) - print(f"✓ Cache saved and index updated") + print("✓ Cache saved and index updated") return True except Exception as e: @@ -641,9 +679,6 @@ def save_to_cache(cache_dir, sofa_file, data_to_cache, interpolation='SH'): # VARIOUS # ----------------------------------- - -import psutil -import os def print_memory_usage(label=""): """Print current memory usage.""" process = psutil.Process(os.getpid()) From 8bbd5cf01f5a8c1fb968a4dcb73a8d474d07e4d6 Mon Sep 17 00:00:00 2001 From: Emanuele Zanoni Date: Tue, 13 Jan 2026 20:36:48 +0100 Subject: [PATCH 6/6] fix(bayesian_listener): revert NPY002 Generator changes - Revert changes to np.random.seed and np.random.multivariate_normal in bayesian_listener.py and utils.py to restore deterministic test behavior. - This is required because pytest test_model_single fails due to non-reproducible random numbers after switching to np.random.Generator. - TODO: properly migrate to Generator API in a way that preserves test determinism. --- bayesian_listener/bayesian_listener.py | 7 ++++--- bayesian_listener/utils.py | 6 +++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/bayesian_listener/bayesian_listener.py b/bayesian_listener/bayesian_listener.py index 95fe4a0..c577abd 100644 --- a/bayesian_listener/bayesian_listener.py +++ b/bayesian_listener/bayesian_listener.py @@ -260,7 +260,8 @@ def infer(self, repetitions = 50, seed = None, prior = 'horizontal'): - rng = np.random.default_rng(seed) + np.random.seed(seed) + # prepare features # use original HRIR if no target is provided @@ -332,7 +333,7 @@ def infer(self, L = np.linalg.cholesky(sigma) # L @ L.T = sigma for t in range(target_num): ts = np.tile(target_feat[t,:], [repetitions, 1]) - xs = ts + rng.normal(size=ts.shape) @ L.T + xs = ts + np.random.normal(size=ts.shape) @ L.T loglik = utils.multiple_logpdfs_vec_input_single_cov( xs,template_feat, logdet, Us).squeeze() logpost = loglik + np.log(prior) @@ -350,7 +351,7 @@ def infer(self, for t in range(target_num): # for ta in range(target_num): # AWGN NOISE - x = rng.multivariate_normal(target_feat[t,:], sigma) + x = np.random.multivariate_normal(target_feat[t,:], sigma) # COMPUTE POSTERIOR # using vectorised solution diff --git a/bayesian_listener/utils.py b/bayesian_listener/utils.py index 8bcbf2f..b8e63bc 100644 --- a/bayesian_listener/utils.py +++ b/bayesian_listener/utils.py @@ -250,7 +250,7 @@ def scatter_von_mises(dirs, sigma_m): return dirs_new def randvmf(kappa, mu, seed = None): - rng = np.random.default_rng(seed) + np.random.seed(seed) assert mu is not None @@ -265,11 +265,11 @@ def randvmf(kappa, mu, seed = None): # Rubinstein 81, p.39, Fisher 87, p.59 kappaS = np.sign(kappa) kappa = abs(kappa) - U = rng.rand() + U = np.random.rand() x = np.log(2. * U * np.sinh(kappa) + np.exp(-kappa)) / kappa x = kappaS * x - psi = 2. * np.pi * rng.rand() + psi = 2. * np.pi * np.random.rand() s_x = np.sqrt(1. - x**2.) y = np.array([np.cos(psi) * s_x, np.sin(psi) * s_x, x])