Conversation
|
Remove python 3.9 ? very old version.
Packaging — ValidationTwo notebooks under
Numbers worth quoting:
The planar kernel is checked for non-regression throughout, including its masked Behaviour changes reviewers should know
Known limitations
Follow-up, not addressed here
|
| ########################################################################### | ||
| def __init__(self, array, nside=None, cell_ids=None, nest=True): | ||
| @classmethod | ||
| def _infer_N0(cls, array): |
There was a problem hiding this comment.
_infer_N0 n'est pas appliquée dans ce fichier. L'inférence de N0 est faite par _infer_N0 du fichier base_dataclass. Cependant, celle-ci renvoie le nombre de pixels (Npix,) au lieu de déterminer plutôt le paramètre nside (dans le cas idéal full-sky).
There was a problem hiding this comment.
Signe « − » sur le troisième terme dans vec_np[:, 2] = np.sqrt(1.0 - vec_np[:, 0] ** 2 + vec_np[:, 1] ** 2).
There was a problem hiding this comment.
Import interne dans une méthode à éviter : torch et numpy sont déjà importés au niveau du module. Remarque identique pour les méthodes _is_varlength_batch (l. 1352), Convol_torch (l. 951) et make_matrix (l. 1144).
There was a problem hiding this comment.
Si, pour une jauge
There was a problem hiding this comment.
Oui, tu as raison. Dans cette ancienne implémentation, si aucun des quatre voisins n’était présent, present[:, col] était entièrement faux. Le traitement de zero_cols faisait alors :
w[0, zero_cols[0]] = present[0, zero_cols[0]].to(w.dtype)
Comme present[0, col] valait également False, cette affectation remettait simplement un poids nul. La somme de la colonne restait donc nulle. Le clamp_min(1e-12) empêchait une division par zéro, mais le résultat de l’interpolation restait une contribution entièrement nulle. Ce n’était donc pas un véritable fallback.
Le code tentait auparavant de remplacer une colonne vide par la colonne correspondant au centre du stencil. Cela traitait le cas courant où le centre possédait au moins un voisin présent dans le domaine local, mais pas le cas où cette colonne centrale était elle-même vide.
Cette logique n’est désormais plus utilisée : l’implémentation locale de SphericalStencil a été supprimée et les convolutions ainsi que la gestion de leur support sont déléguées à healpix-analyse. Il n’y a donc plus, dans STL, cette normalisation locale des poids ni ce cas silencieux où une colonne vide produit une sortie nulle.
| src = ref_cols[empty_cols] | ||
|
|
||
| # copie idx/w de la colonne 'centre' | ||
| # copy idx/w from the 'centre' column |
There was a problem hiding this comment.
Prévoir un plan de secours si le centre du stencil local n'a lui aussi aucun voisin présent ?
There was a problem hiding this comment.
Oui, un second niveau de fallback aurait été nécessaire dans cette implémentation.
La copie de la colonne centrale supposait implicitement que le centre du stencil avait au moins un voisin présent dans ids_sorted. Si cette hypothèse était fausse, la colonne copiée restait vide et la normalisation suivante produisait uniquement des poids nuls. Le calcul restait numériquement fini grâce au clamp_min, mais le résultat n’avait pas de signification correcte : il revenait à remplacer silencieusement la contribution par zéro.
Si nous avions conservé cette implémentation, le comportement robuste aurait été :
- tenter d’utiliser les voisins du point demandé ;
- en leur absence, tenter le centre du stencil ;
- si le pixel cible appartient au domaine local, imposer explicitement une interpolation identité sur ce pixel avec un poids égal à un ;
- si même le pixel cible n’appartient pas au domaine, signaler explicitement que la sortie n’est pas définie, plutôt que de produire silencieusement zéro.
Toutefois, ce code de stencil local faisait double emploi avec healpix-analyse. En supprimant SphericalStencil de STL et utilisant directement les opérateurs de healpix-analyse, qui portent maintenant la responsabilité de la géométrie, de l’interpolation et du support partiel. Ce cas limite n’a donc plus à être corrigé dans STL.
| Keep only the frequencies inside the Nyquist radius of the pixel grid. | ||
| """ | ||
| from STL_main.Synthesis import apply_nyquist_filter | ||
|
|
There was a problem hiding this comment.
Mettre cet import en haut du fichier
| data.array = torch.abs(data.array) | ||
| npix_full = 12 * self.N0[0] ** 2 | ||
| if array.shape[-1] != npix_full: | ||
| return array |
There was a problem hiding this comment.
Mettre un petit warning pour prévenir, lors de la synthèse, que le filtrage bandlimit n'a pas pu être effectué à cause du masque.
| [project] | ||
| name = "STL" | ||
| dynamic = ["version"] | ||
| requires-python = ">=3.10" |
There was a problem hiding this comment.
À l’origine de l’échec du job test(3.9), faut-il retirer ce test ou abaisser l’exigence de version de Python du projet ?
There was a problem hiding this comment.
Il faut retirer Python 3.9 de la matrice de CI plutôt que d’abaisser l’exigence de version du projet.
STL déclare maintenant :
requires-python = ">=3.10"
et healpix-analyse, qui est devenu une dépendance directe de STL pour le noyau HEALPix, déclare également requires-python = ">=3.10". Le job Python 3.9 échoue donc avant même de tester réellement le code : pip refuse à juste titre d’installer un projet dont la version minimale annoncée est Python 3.10.
Conserver le job 3.9 tout en déclarant Python 3.10 comme version minimale est contradictoire. Abaisser seulement la contrainte de STL ne suffirait pas non plus, puisque la dépendance healpix-analyse resterait incompatible avec Python 3.9. Rétablir une prise en charge réelle de Python 3.9 demanderait d’abaisser la contrainte dans les deux projets, puis de vérifier toutes les dépendances et l’ensemble des tests sous cette version. Cela dépasse le périmètre de cette PR.
La correction cohérente consiste donc à supprimer 3.9 de .github/workflows/ci.yaml et à conserver les jobs à partir de Python 3.10. La version minimale déclarée par le paquet et les versions testées par la CI seront ainsi alignées.
Le 3.9 était lié a de vieux driver GPU NVIDIA, je pense que c'est fini.
| # Reference is the head of the sibling checkout, installed in editable mode: | ||
| # pip install -e ../healpix-analyse | ||
| "healpix-analyse", | ||
| ] |
There was a problem hiding this comment.
Redondance des dépendances avec requirements.txt : tout regrouper dans pyproject.toml.
| "healpix-resample", | ||
| "xarray", | ||
| "zarr", | ||
| ] |
There was a problem hiding this comment.
Ajout des dépendances optionnelles provenant de requirements.txt pour les tests et la documentation :
[project.optional-dependencies]
dev = [
"pytest",
"pre-commit",
]
docs = [
"sphinx",
"sphinx-rtd-theme",
"sphinx-autodoc-typehints",
]
| import torch.nn.functional as F | ||
| from healpix_analyse.convol import HealPixConv | ||
| from healpix_analyse.down import HealPixDown | ||
| from healpix_analyse.healpix_sht import HEALPixSHT |
There was a problem hiding this comment.
Pourquoi laisser SphericalStencil dans le dossier source si on utilise healpix-analyse ? Et healpix-analyse est-elle basée sur la classe SphericalStencil pour ses convolutions, downgrading, etc. ?
There was a problem hiding this comment.
Tu as raison : une fois healpix-analyse utilisé comme source de la géométrie sphérique, il n’y avait plus de raison de conserver une seconde implémentation de SphericalStencil dans STL.
SphericalStencil correspondait à une ancienne implémentation interne qui gérait notamment la recherche des voisins, les poids d’interpolation, les domaines partiels et les opérations de montée ou de descente en résolution. La conserver en parallèle de healpix-analyse créait deux implémentations de la même logique, avec le risque qu’elles divergent ou que des corrections soient appliquées à l’une mais pas à l’autre. Les deux questions précédentes sur les colonnes sans voisin illustrent précisément ce problème.
Le fichier STL_main/SphericalStencil.py a donc été supprimé. Le noyau HEALPix de STL importe maintenant directement :
from healpix_analyse.convol import HealPixConv
from healpix_analyse.down import HealPixDown
from healpix_analyse.healpix_sht import HEALPixSHT
healpix-analyse n’est pas basé sur la classe SphericalStencil de STL. C’est une implémentation indépendante :
- HealPixConv construit et applique sa propre géométrie de convolution, ses poids d’interpolation et ses jauges ;
- HealPixDown effectue la réduction de résolution avec son propre opérateur ;
- HEALPixSHT fournit les transformations harmoniques utilisées pour le filtrage en bande et les spectres angulaires ;
- les domaines locaux sont décrits directement par cell_ids.
STL conserve uniquement la couche d’adaptation nécessaire au calcul des statistiques de scattering : création et mise en cache des opérateurs, préparation des kernels complexes, suivi de conv_history, propagation des masques et intégration avec ST_Operator.
La dépendance est donc dans un seul sens : STL utilise les opérateurs publics de healpix-analyse; healpix-analyse ne dépend pas de SphericalStencil ni de STL. Cela évite la duplication et permet de centraliser les futures corrections géométriques dans healpix-analyse.
| raise Exception( | ||
| f"Data should be a STL_Healpix_Kernel_Torch instance, got {type(data)}" | ||
| ) | ||
| if self.DT != data.DT: |
There was a problem hiding this comment.
Rendre DT attribut de classe de WaveletOperatorHealpixKernel_Torch (faire de même pour les autres opérateurs wavelet). Cela évitera d’instancier des opérateurs wavelet avec des DT divergents de ceux des classes data correspondantes, qui sont déjà des attributs de classe.
class WaveletOperatorHealpixKernel_Torch:
DT: ClassVar[str] = "HealpixKernel_Torch"Le second if devient alors redondant.
|
|
||
| ########################################################################### | ||
| def get_L(self): | ||
| return self.L | ||
|
|
| def _smooth_with_nan(self, data, inplace: bool = True): | ||
| """ | ||
| NaN-aware smoothing: the map and the validity mask are both smoothed, | ||
| and the result is normalized by the sum of the valid weights. |
There was a problem hiding this comment.
Indiquer en docstring :
"""The last dimension of `data.array` corresponds to `K`, the total number of pixels in the full-sky grid at this resolution. Pixels with no data are represented by `NaN` values."""car, en général, la dernière dimension d’une instance data est de taille N_pix, avec N_pix <= K en présence d’un masque.
There was a problem hiding this comment.
Attention, STL Healpix peut fonctionner avec une liste de cell_ids healpix pour travailler en local. Il ne faut pas casser cette propriété.
| f"Npix pixels), got {tuple(m.shape)}." | ||
| ) | ||
|
|
||
| invalid = m.to(torch.bool) if m.dtype == torch.bool else (m != 0) |
There was a problem hiding this comment.
invalid = m if m.dtype == torch.bool else (m != 0)
| ) | ||
|
|
||
| self.bin_windows = windows # [n_bins, lmax+1] | ||
| self.bin_weights = windows * (2.0 * ell + 1.0)[None, :] |
There was a problem hiding this comment.
Permet d'obtenir un estimateur de moyenne pondérée
| """ | ||
| Build the [n_bins, lmax+1] window matrix and its (2l+1) weighting. | ||
| """ | ||
| l_min = max(1.0, float(2**self.Jmin)) |
There was a problem hiding this comment.
Pourrait être intéressant de discuter des bornes du spectre de puissance et de voir s’il serait pertinent de les faire dépendre de l’échelle
| # the monopole carries no information and l < l_min is excluded | ||
| windows = torch.where( | ||
| (ell >= l_min)[None, :], windows, torch.zeros_like(windows) | ||
| ) |
There was a problem hiding this comment.
Pourrait on obtenir un gain temps/mémoire intéressant en ne calculant pas les modes dont le multipôle est inférieur à
| x_flat = x.reshape(B, K) | ||
| full_sky = bool(data.pbc) | ||
| if use_band_maps is None: | ||
| use_band_maps = not full_sky |
There was a problem hiding this comment.
Avertir de la présence d'artefacts si on a une carte partiel et que l'on sélectionne tout de même l'option use_band_maps=False
if not use_band_maps and not full_sky:
warnings.warn(
"`use_band_maps=False` with partial-sky data: "
"missing pixels are zero-filled before computing the cross spectrum "
"using `anafast` or `alm`, which may introduce artifacts."
)
| data.nside = data.N0[0] // (2**dg_out) | ||
| ########################################################################### | ||
| def _apply_anafast(self, maps_full, pairs): | ||
| """Full-sky estimator going through HEALPixSHT.anafast pair by pair.""" |
There was a problem hiding this comment.
Ne risque-t-on pas de recalculer plusieurs fois les coefficients HEALPixSHT.anafast ?
2 participants
What this does
Brings
STL_Healpix_Kernel_Torchup to the same interface asSTL_2D_Kernel_Torch, so that the data-type independent machinery(
ST_Operator,ST_Statistics,Synthesis) runs unchanged on the sphere.It also removes the
foscatdependency:SphericalStencilimported it at moduleload, which made
import STL_main.STL_Healpix_Kernel_Torchfail outright on anyenvironment without that package. The spherical geometry now comes from
healpix-analyse.Why the HEALPix kernel had diverged
Not in the spherical convolution, but in the division of responsibilities: the
HEALPix version kept the statistics (
mean,cov) on the data class, whilethe planar kernel had moved them onto the wavelet operator — the only object
that knows the mask, the convolution history and the borders.
Changes
Data-type independence
Base_DataClass:NDIM_PIX(2 planar, 1 HEALPix) and an overridable_infer_N0;new_like()to build a companion field on the same geometry;apply_bandlimit()as a hook for the band-limiting of an initial random field.N0is no longer overwritten when it is passed explicitly.ST_Operator.apply: six slicings assumed two pixel axes(
data_l1m[j3][:, ch, :, :, :]). The trailing:were redundant; droppingthem makes the code work for any number of pixel axes, with no change in 2D.
torch_backend.maskmean: accepts any set of trailing dimensions, not only(-2, -1).ST_Statistics: carriesdata_example(an empty copy of the input) andpix_shape, which is notN0for every data type.HEALPix kernel
Base_DataClass, withpbc/dg/N0=(nside,)/conv_history/cell_ids/nest; gainsdivide,get_ST_op,get_CS_op.mean,square_mean,cov,standardize,unstandardize,_compute_and_store_cross_cov,j_to_dg,mask_full_res, and adownsamplewith the planar signature.healpix_analyse.convol.HealPixConv: the complex kernel iscarried as two output channels, the
Lorientations are theLgauges, so onecall returns the complex answer. Anti-aliased decimation via
healpix_analyse.down.HealPixDown, one level at a time.canonical order takes the package's fast path.
Angular power spectrum (
CS_operator_Healpix_Torch)n_bins,bin_centers,apply(...) -> [Nb, Nc, Nc, n_bins],plot_cross_spectrum), but the estimatoris
C_ell, built onhealpix_analyse.healpix_sht.HEALPixSHT.map2almonce per channel, then every requested pair —Nctransforms instead of
Nc². A reference route throughanafastpair by pairis kept behind
cross_spectrum_method="anafast".cell_ids, band-filtered withalm2map, and the cross-covariance is taken over the observed pixels,C_b = 4π ⟨f_b·g⟩_obs / Σ_l (2l+1) W_b(l).(2l+1)times a log-Gaussianwindow — the spherical transposition of
_build_log_gaussian_bin_masks.Masked data
precomputed once, the spherical counterpart of
_build_reweighting_maps_and_scattering_layer_masks. The mask is eroded by akernel of ones over the
Lgauges, i.e. the exact union of the supports used;HealPixDownapplied to the mask gives the local invalid fractionf, and acoarse pixel is either declared invalid past
downsample_nan_weight_thresholdor rescaled by1/(1-f).mask_full_res=Falseopts out.cov: for S4 both operands sit at the same depth but ondifferent scales (
|I*ψ_j1|*ψ_j3vs|I*ψ_j2|*ψ_j3), so neither mask containsthe other and the union is required. Only the first was used.
SynthesisNDIM_PIX; companion fields built from a prototypethrough
new_likerather than from the class alone; the Nyquist prefilterbecomes
apply_bandlimiton the data class (Nyquist disc in the plane,spherical harmonic round trip on the sphere).
ScatteringMatchModelacceptseither a class (legacy) or a prototype, so the planar path is unchanged.
SphericalStencilimport foscat.scat_covremoved.Down/Upreimplemented in healpy +torch:
Downreproduceshp.ud_gradeexactly,Upreproduceshp.get_interp_valto 3e-16, both differentiable and complex-capable. Alsofixes
self.cell_ids→self.cell_ids_default, which made thecell_ids=Nonebranch of both raiseAttributeError.scat_opconstructor argument