foscat.scat_cov =============== .. py:module:: foscat.scat_cov Attributes ---------- .. autoapisummary:: foscat.scat_cov.tf_defined foscat.scat_cov.tf_function foscat.scat_cov.testwarn Classes ------- .. autoapisummary:: foscat.scat_cov.scat_cov foscat.scat_cov.funct Functions --------- .. autoapisummary:: foscat.scat_cov.read Module Contents --------------- .. py:data:: tf_defined .. py:data:: tf_function .. py:function:: read(filename) .. py:data:: testwarn :value: 0 .. py:class:: scat_cov(s0, s2, s3, s4, s1=None, s3p=None, backend=None, use_1D=False, return_data=False) .. py:attribute:: S0 .. py:attribute:: S2 .. py:attribute:: S3 .. py:attribute:: S4 .. py:attribute:: S1 :value: None .. py:attribute:: S3P :value: None .. py:attribute:: backend :value: None .. py:attribute:: idx1 :value: None .. py:attribute:: idx2 :value: None .. py:attribute:: use_1D :value: False .. py:method:: numpy() .. py:method:: constant() .. py:method:: conv2complex(val) .. py:method:: flatten() .. py:method:: flattenMask() .. py:method:: get_S0() .. py:method:: get_S1() .. py:method:: get_S2() .. py:method:: reset_S2() .. py:method:: get_S3() .. py:method:: get_S3P() .. py:method:: get_S4() .. py:method:: get_j_idx() .. py:method:: get_js4_idx() .. py:method:: relu() .. py:method:: domult(x, y) .. py:method:: dodiv(x, y) .. py:method:: domin(x, y) .. py:method:: doadd(x, y) .. py:method:: interp(nscale, extend=True, constant=False) .. py:method:: plot(name=None, hold=True, color='blue', lw=1, legend=True, norm=False) .. py:method:: get_np(x) .. py:method:: save(filename) .. py:method:: read(filename) .. py:method:: std() .. py:method:: mean() .. py:method:: initdx(norient) .. py:method:: sqrt() .. py:method:: L1() .. py:method:: square_comp() .. py:method:: iso_mean(repeat=False) .. py:method:: fft_ang(nharm=1, imaginary=False) Project orientation axes onto the first Fourier harmonics. This is a softer alternative to :meth:`iso_mean`. Instead of collapsing each orientation axis to a single mean value, it keeps the first ``nharm`` harmonics of the discrete Fourier transform along each orientation axis L. This preserves the *amplitude* of the angular variation while reducing the number of descriptors from L to ``nout = 1 + 2*nharm`` (with ``imaginary=True``) per orientation axis. .. rubric:: Why ``imaginary=True`` is the recommended mode With ``imaginary=False`` the projection basis is ``{1, cos(2π·l/L), cos(4π·l/L), …}``. A field whose dominant orientation sits at exactly 90° (i.e. at the zero-crossing of the cosine) would give a near-zero first-harmonic coefficient even though the field is strongly anisotropic. With ``imaginary=True`` the basis is ``{1, cos(2π·l/L), sin(2π·l/L), cos(4π·l/L), sin(4π·l/L), …}``. The amplitude of the first harmonic, .. math:: A_1 = \sqrt{c_1^2 + s_1^2} is **rotation-invariant**: it is the same regardless of the absolute orientation of the field. Use ``imaginary=True`` whenever you want a description that does not depend on the image orientation. .. rubric:: Reduction per statistic ============ ==================== ========================== Statistic Input shape Output shape (nharm=1) ============ ==================== ========================== S1, S2 ``[…, L]`` ``[…, nout]`` S3, S3P ``[…, L, L]`` ``[…, L, nout]`` S4 ``[…, L, L, L]`` ``[…, L, L, nout]`` ============ ==================== ========================== For S1/S2 the projection applies the Fourier basis directly on the single orientation axis. For S3/S3P and S4 the projection is **not** a tensor product. Instead, the statistics are first reindexed by the relative-orientation differences (exactly as in :meth:`iso_mean`), then the absolute orientation axis :math:`l_1` is projected onto the Fourier basis: .. math:: \text{S3\_out}[\Delta l,\, k] = \sum_{l_1} \phi_k(l_1)\, S3[l_1,\,(l_1+\Delta l)\bmod L] .. math:: \text{S4\_out}[\Delta l_{12},\, \Delta l_{13},\, k] = \sum_{l_1} \phi_k(l_1)\, S4[l_1,\,(l_1+\Delta l_{12})\bmod L,\,(l_1+\Delta l_{13})\bmod L] where :math:`\phi_0(l)=1`, :math:`\phi_1(l)=\cos(2\pi l/L)`, :math:`\phi_2(l)=\sin(2\pi l/L)`. This preserves the relative-orientation axes (same as iso_mean) while also capturing how strongly the statistics vary as the whole frame rotates. The :math:`k=0` component is identical to the :meth:`iso_mean` result. :Parameters: * **nharm** (:py:class:`int`, *optional*) -- Number of harmonics to keep beyond the DC term. ``nharm=1`` (default) keeps the mean and the first angular harmonic. * **imaginary** (:py:class:`bool`, *optional*) -- If ``False`` (default), keep only the cosine components ``{1, cos, cos(2·), …}`` — ``nout = 1 + nharm``. If ``True``, keep both cosine and sine components ``{1, cos, sin, cos(2·), sin(2·), …}`` — ``nout = 1 + 2·nharm``. **Recommended: ``True``** to obtain rotation-invariant amplitudes. :returns: :py:class:`scat_cov` -- A new statistics object with orientation axes compressed from L to ``nout``. .. admonition:: Examples >>> stat = scat_op.eval(image) >>> stat_fft = stat.fft_ang(nharm=1, imaginary=True) >>> # S2 shape: (..., L) -> (..., 3) [DC, cos, sin] >>> # S3 shape: (..., L, L) -> (..., L, 3) [Δl axis kept, l1 compressed] >>> # S4 shape: (..., L, L, L) -> (..., L, L, 3) >>> >>> # Rotation-invariant first-harmonic amplitude for S2: >>> import numpy as np >>> A1_S2 = np.sqrt(stat_fft.S2[..., 1]**2 + stat_fft.S2[..., 2]**2) >>> >>> # Angular modulation of S3 at a given relative angle Δl: >>> A1_S3_delta0 = np.sqrt(stat_fft.S3[..., 0, 1]**2 + stat_fft.S3[..., 0, 2]**2) .. py:method:: fft_ang_sigma(nharm=1, imaginary=False) Propagate per-coefficient standard deviations through ``fft_ang``. When the input object contains standard deviations ``sigma_j`` (e.g. the variance output of :meth:`eval` with ``calc_var=True``), applying :meth:`fft_ang` directly would be **wrong** because linear projection of sigma gives near-zero values for harmonic components (cosine/sine sum to zero over uniform sigma values). The correct formula for a linear map :math:`y_k = \sum_j A_{jk}\,x_j` with independent inputs of variance :math:`\sigma_j^2` is: .. math:: \sigma_{y_k} = \sqrt{\sum_j A_{jk}^2 \, \sigma_j^2} This method applies that formula element-wise using the squared projection matrices, ensuring that the output sigma is always positive and physically meaningful for use as the ``sigma`` argument of :meth:`reduce_distance`. :Parameters: * **nharm** (:py:class:`int`, *optional*) -- Passed through to the underlying projection matrices (same as in :meth:`fft_ang`). * **imaginary** (:py:class:`bool`, *optional*) -- Passed through (same as in :meth:`fft_ang`). :returns: :py:class:`scat_cov` -- A new statistics object with the same shape as ``fft_ang`` output but containing propagated standard deviations. .. py:method:: iso_std(repeat=False) .. py:method:: get_nscale() .. py:method:: get_norient() .. py:method:: add_data_from_log_slope(y, n, ds=3) .. py:method:: add_data_from_slope(y, n, ds=3) .. py:method:: up_grade(nscale, ds=3) .. py:class:: funct(NORIENT=4, LAMBDA=1.2, KERNELSZ=3, slope=1.0, all_type='float32', nstep_max=20, padding='SAME', gpupos=0, mask_thres=None, mask_norm=False, isMPI=False, TEMPLATE_PATH=None, BACKEND='torch', use_2D=False, use_1D=False, return_data=False, DODIV=False, use_median=False, InitWave=None, silent=True, mpi_size=1, mpi_rank=0) Bases: :py:obj:`foscat.FoCUS.FoCUS` .. py:method:: fill(im, nullval=hp.UNSEEN) .. py:method:: moments(list_scat) .. py:method:: stat_cfft(im, image2=None, upscale=False, smooth_scale=0, spin=0) .. py:method:: stat_cfft_cell_ids(im, image2=None, upscale=False, smooth_scale=0, spin=0, cell_ids=None) Full-sky stat_cfft + per-scale slicing by cell_ids. If cell_ids is None, identical to stat_cfft. If provided, cell_ids must be NEST indices at the input resolution. .. py:method:: div_norm(complex_value, float_value) .. py:method:: eval(image1, image2=None, mask=None, norm=None, calc_var=False, cmat=None, cmat2=None, Jmax=None, out_nside=None, edge=True, nside=None, cell_ids=None, spin=0) Calculates the scattering correlations for a batch of images. Mean are done over pixels. mean of modulus: S1 = <|I * Psi_j3|> Normalization : take the log power spectrum: S2 = <|I * Psi_j3|^2> Normalization : take the log orig. x modulus: S3 = < (I * Psi)_j3 x (|I * Psi_j2| * Psi_j3)^* > Normalization : divide by (S2_j2 * S2_j3)^0.5 modulus x modulus: S4 = <(|I * psi1| * psi3)(|I * psi2| * psi3)^*> Normalization : divide by (S2_j1 * S2_j2)^0.5 :Parameters: * **image1** (:py:class:`tensor`) -- Image on which we compute the scattering coefficients [Nbatch, Npix, 1, 1] * **image2** (:py:class:`tensor`) -- Second image. If not None, we compute cross-scattering covariance coefficients. * **mask** * **norm** (:py:obj:`None` or :py:class:`str`) -- If None no normalization is applied, if 'auto' normalize by the reference S2, if 'self' normalize by the current S2. * **spin** (:py:class:`Integer`) -- If different from 0 compute spinned data (U,V to Divergence/Rotational spin==1) or (Q,U to E,B spin=2). This implies that the input data is 2*12*nside^2. :returns: :py:class:`S1`, :py:class:`S2`, :py:class:`S3`, :py:class:`S4 normalized` .. py:method:: clean_norm() .. py:method:: computer_filter(M, N, J, L) This function is strongly inspire by the package https://github.com/SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel. .. py:method:: cut_high_k_off(data_f, dx, dy) This function is strongly inspire by the package https://github.com/SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel. .. py:method:: get_dxdy(j, M, N) This function is strongly inspire by the package https://github.com/SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel. .. py:method:: get_edge_masks(M, N, J, d0=1, in_mask=None, edge_dx=None, edge_dy=None) This function is strongly inspire by the package https://github.com/SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel. .. py:method:: scattering_cov(data, data2=None, Jmax=None, if_large_batch=False, S4_criteria=None, use_ref=False, normalization='S2', edge=False, in_mask=None, pseudo_coef=1, get_variance=False, ref_sigma=None, iso_ang=False, return_table=False, fft_ang=False, fft_nharm=1, fft_imaginary=True) Calculates the scattering correlations for a batch of images, including: This function is strongly inspire by the package https://github.com/SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel. orig. x orig.: P00 = <(I * psi)(I * psi)*> = L2(I * psi)^2 orig. x modulus: C01 = <(I * psi2)(|I * psi1| * psi2)*> / factor when normalization == 'P00', factor = L2(I * psi2) * L2(I * psi1) when normalization == 'P11', factor = L2(I * psi2) * L2(|I * psi1| * psi2) modulus x modulus: C11_pre_norm = <(|I * psi1| * psi3)(|I * psi2| * psi3)> C11 = C11_pre_norm / factor when normalization == 'P00', factor = L2(I * psi1) * L2(I * psi2) when normalization == 'P11', factor = L2(|I * psi1| * psi3) * L2(|I * psi2| * psi3) modulus x modulus (auto): P11 = <(|I * psi1| * psi2)(|I * psi1| * psi2)*> :Parameters: * **data** (:py:class:`numpy array` or :py:class:`torch tensor`) -- image set, with size [N_image, x-sidelength, y-sidelength] * **if_large_batch** (:py:class:`Bool (=False)`) -- It is recommended to use "False" unless one meets a memory issue * **C11_criteria** (:py:class:`str` or :py:class:`None (=None)`) -- Only C11 coefficients that satisfy this criteria will be computed. Any expressions of j1, j2, and j3 that can be evaluated as a Bool is accepted.The default "None" corresponds to "j1 <= j2 <= j3". * **use_ref** (:py:class:`Bool (=False)`) -- When normalizing, whether or not to use the normalization factor computed from a reference field. For just computing the statistics, the default is False. However, for synthesis, set it to "True" will stablize the optimization process. * **normalization** (str ``'P00'`` or ``'P11'`` (=``'P00'``:py:class:`)`) -- Whether 'P00' or 'P11' is used as the normalization factor for C01 and C11. * **remove_edge** (:py:class:`Bool (=False)`) -- If true, the edge region with a width of rougly the size of the largest wavelet involved is excluded when taking the global average to obtain the scattering coefficients. :returns: * **'P00'** (:py:class:`torch tensor with size [N_image`, ``J``, :py:class:`L] (# image`, :py:class:`j1`, :py:class:`l1)`) -- the power in each wavelet bands (the orig. x orig. term) * **'S1'** (:py:class:`torch tensor with size [N_image`, ``J``, :py:class:`L] (# image`, :py:class:`j1`, :py:class:`l1)`) -- the 1st-order scattering coefficients, i.e., the mean of wavelet modulus fields * **'C01'** (:py:class:`torch tensor with size [N_image`, ``J``, ``J``, :py:class:`L`, :py:class:`L] (# image`, :py:class:`j1`, :py:class:`j2`, :py:class:`l1`, :py:class:`l2)`) -- the orig. x modulus terms. Elements with j1 < j2 are all set to np.nan and not computed. * **'C11'** (:py:class:`torch tensor with size [N_image`, ``J``, ``J``, ``J``, :py:class:`L`, :py:class:`L`, :py:class:`L] (# image`, :py:class:`j1`, :py:class:`j2`, :py:class:`j3`, :py:class:`l1`, :py:class:`l2`, :py:class:`l3)`) -- the modulus x modulus terms. Elements not satisfying j1 <= j2 <= j3 and the conditions defined in 'C11_criteria' are all set to np.nan and not computed. * **'C11_pre_norm' and 'C11_pre_norm_iso'** (:py:class:`pre-normalized modulus x modulus terms.`) * **'P11'** (:py:class:`torch tensor with size [N_image`, ``J``, ``J``, :py:class:`L`, :py:class:`L] (# image`, :py:class:`j1`, :py:class:`j2`, :py:class:`l1`, :py:class:`l2)`) -- the modulus x modulus terms with the two wavelets within modulus the same. Elements not following j1 <= j3 are set to np.nan and not computed. * **'P11_iso'** (:py:class:`torch tensor with size [N_image`, ``J``, ``J``, :py:class:`L] (# image`, :py:class:`j1`, :py:class:`j2`, :py:class:`l2-l1)`) -- 'P11' averaged over l1 while keeping l2-l1 constant. .. py:method:: purge_edge_mask() .. py:method:: to_gaussian(x, in_mask=None) .. py:method:: from_gaussian(x) .. py:method:: square(x) .. py:method:: sqrt(x) .. py:method:: reduce_mean(x) .. py:method:: reduce_mean_batch(x) .. py:method:: reduce_sum_batch(x) .. py:method:: reduce_distance(x, y, sigma=None) .. py:method:: reduce_sum(x) .. py:method:: ldiff(sig, x) .. py:method:: log(x) .. py:method:: eval_comp_fast(image1, image2=None, mask=None, norm=None, cmat=None, cmat2=None) .. py:method:: eval_fast(image1, image2=None, mask=None, norm=None, cmat=None, cmat2=None) .. py:method:: calc_matrix_orientation(noise_map, image2=None) .. py:method:: synthesis(image_target, reference=None, nstep=4, seed=1234, Jmax=None, edge=False, to_gaussian=False, use_variance=True, synthesised_N=1, input_image=None, grd_mask=None, in_mask=None, iso_ang=False, fft_ang=False, fft_nharm=1, fft_imaginary=True, EVAL_FREQUENCY=100, NUM_EPOCHS=300, scat_cov_method='eval', n_up=0) Synthesise a new field whose scattering-covariance statistics match those of ``image_target``. This is the main high-level entry point for texture synthesis, denoising, inpainting, and component separation with FOSCAT. Internally it runs L-BFGS-B (``scipy.optimize.fmin_l_bfgs_b``) on the scattering-covariance loss, with an optional coarse-to-fine multi-resolution schedule controlled by ``nstep``. The optimisation minimises: .. math:: u^* = \arg\min_u \mathcal{L}(u), \qquad \mathcal{L}(u) = \sum_k \frac{(\Phi(u)_k - \Phi(d)_k)^2}{\sigma_k^2} where :math:`\Phi` is the scattering-covariance operator, :math:`d` is ``image_target``, and :math:`\sigma_k^2` is the per-coefficient variance of :math:`\Phi(d)` (used only when ``use_variance=True``). :Parameters: * **image_target** (:py:class:`array-like`) -- Reference field whose statistics are to be reproduced. - **HEALPix:** shape ``(npix,)`` or ``(B, npix)`` with ``npix = 12 * nside**2``. - **2D:** shape ``(H, W)`` or ``(B, H, W)``. - **1D:** shape ``(N,)`` or ``(B, N)``. Normalise the input before calling (zero mean, unit variance is standard):: xnorm = (image - image.mean()) / image.std() * **reference** (:py:class:`array-like`, *optional*) -- Second reference field for **cross-covariance** synthesis. When provided the loss matches the cross-statistics :math:`\Phi_\times(u, d_2)` against the target :math:`\Phi_\times(d_1, d_2)`, where ``d1 = image_target`` and ``d2 = reference``. Both fields must have the same shape. Useful for component separation (e.g. CMB × dust template). * **nstep** (:py:class:`int`, *optional*) -- Number of resolution levels in the coarse-to-fine schedule. Default is ``4``. The algorithm downsamples ``image_target`` by successive factors of 2 to build a resolution pyramid, then optimises from coarsest to finest, using each result as the warm start for the next level. This dramatically improves convergence for large maps. +---------+----------------------------------------------+ | nstep | Resolutions visited (2D, 256 × 256 target) | +=========+==============================================+ | 1 | 256 × 256 only | +---------+----------------------------------------------+ | 2 | 128 × 128 → 256 × 256 | +---------+----------------------------------------------+ | 3 | 64 × 64 → 128 × 128 → 256 × 256 | +---------+----------------------------------------------+ | 4 | 32 → 64 → 128 → 256 | +---------+----------------------------------------------+ Capped automatically when ``nstep > jmax - 1`` (map too small for that many downsampling steps). * **seed** (:py:class:`int`, *optional*) -- Random seed for the Gaussian white-noise initial condition used at the coarsest resolution level. Default is ``1234``. Has no effect when ``input_image`` is provided. Change to generate independent realisations:: results = [scat_op.synthesis(xnorm, seed=s) for s in range(8)] * **Jmax** (:py:class:`int` or :py:obj:`None`, *optional*) -- Maximum wavelet scale index included in the loss. ``None`` (default) uses all scales available for the map size. Decremented by 1 at each coarser resolution level during the multi-resolution schedule. Reduce to constrain only small-scale statistics:: result = scat_op.synthesis(xnorm, Jmax=4) # ignore large scales * **edge** (:py:class:`bool`, *optional*) -- Enable edge-aware boundary handling for non-periodic maps. Default is ``False`` (periodic/full-sphere boundaries). Set to ``True`` for rectangular images or partial-sky patches. Activated automatically when ``in_mask`` is provided. * **to_gaussian** (:py:class:`bool`, *optional*) -- Gaussianise ``image_target`` before computing target statistics, then invert at the end. Default is ``False``. Useful for highly non-Gaussian fields (log-normal distributions, etc.) to decouple the histogram constraint from the scattering statistics. * **use_variance** (:py:class:`bool`, *optional*) -- Weight each loss term by the inverse variance of the corresponding coefficient in ``image_target``. Default is ``True``. - ``True``: loss is scale-invariant — large-scale and small-scale coefficients contribute equally regardless of amplitude. - ``False``: uniform weighting, dominated by the largest-amplitude statistics (usually coarse scales). * **synthesised_N** (:py:class:`int`, *optional*) -- Number of independent synthetic maps to produce simultaneously. Default is ``1``. All ``synthesised_N`` maps share the same loss, which amortises the per-iteration GPU cost. The output gains a leading batch dimension:: result = scat_op.synthesis(xnorm, synthesised_N=4) # result.shape == (4, H, W) for 2D * **input_image** (:py:class:`array-like` or :py:obj:`None`, *optional*) -- Warm-start initial condition. When ``None`` (default), starts from Gaussian white noise. When provided, ``input_image`` is downsampled to each resolution level and used as the initial guess at the coarsest level. Useful for iterative refinement or for injecting a prior:: result_v1 = scat_op.synthesis(xnorm, NUM_EPOCHS=100) result_v2 = scat_op.synthesis(xnorm, NUM_EPOCHS=500, input_image=result_v1) * **grd_mask** (:py:class:`array-like` or :py:obj:`None`, *optional*) -- Binary mask selecting which pixels the optimiser is allowed to update. Pixels where ``grd_mask = 0`` are frozen to their values in ``image_target``; only pixels where ``grd_mask = 1`` move. Same shape as ``image_target``. Downsampled at each resolution level. Useful for **inpainting** (reconstruct missing regions while keeping observed pixels fixed):: grd_mask = np.zeros_like(xnorm) grd_mask[missing_pixels] = 1.0 result = scat_op.synthesis(xnorm, grd_mask=grd_mask) * **in_mask** (:py:class:`array-like` or :py:obj:`None`, *optional*) -- Binary mask marking **invalid pixels in the input data** that should not contribute to the reference statistics. Same shape as ``image_target``. Pixels with ``in_mask = 0`` are excluded from the scattering-covariance computation at every scale level (the mask is downsampled along with the map). Setting ``in_mask`` also enables ``edge=True`` internally. Useful for partial-sky CMB maps or images with missing regions:: mask = np.ones_like(xnorm) mask[bad_pixels] = 0.0 result = scat_op.synthesis(xnorm, in_mask=mask) * **iso_ang** (:py:class:`bool`, *optional*) -- Use isotropically averaged statistics in the loss. Default is ``False``. - ``False``: the full oriented S1–S4 statistics are used. The synthesised field can reproduce anisotropic structures (filaments, oriented textures). - ``True``: statistics are collapsed to their rotationally invariant content via :meth:`iso_mean` before computing the loss. Reduces the number of constraints and accelerates convergence for statistically isotropic fields (e.g. CMB). Angular-reduction summary with ``iso_ang=True``: ============ ==================== ======================= Statistic Input shape Output shape ============ ==================== ======================= S1, S2 ``[…, L]`` ``[…]`` (mean over L) S3, S3P ``[…, L, L]`` ``[…, L]`` (by Δl = l₂−l₁) S4 ``[…, L, L, L]`` ``[…, L, L]`` ============ ==================== ======================= .. note:: ``iso_ang`` and ``fft_ang`` should not be used together. ``iso_ang`` is the harder reduction (mean only); ``fft_ang`` is softer and keeps angular variation. * **fft_ang** (:py:class:`bool`, *optional*) -- Use Fourier-compressed angular statistics in the loss. Default is ``False``. Softer alternative to ``iso_ang``: instead of collapsing each orientation axis to a single mean, keeps the first ``fft_nharm`` Fourier harmonics along each axis. With ``fft_nharm=1, fft_imaginary=True`` (defaults), each orientation axis L is projected to 3 coefficients: - index 0 — DC (mean, same as ``iso_ang``); - index 1 — cosine of first harmonic; - index 2 — sine of first harmonic. The amplitude :math:`A_1 = \sqrt{c_1^2 + s_1^2}` is **rotation-invariant**: it does not depend on the absolute orientation of the image. Shape reduction: ============ ==================== ========================== Statistic Input shape Output shape (nharm=1) ============ ==================== ========================== S1, S2 ``[…, L]`` ``[…, 3]`` S3, S3P ``[…, L, L]`` ``[…, 3, 3]`` S4 ``[…, L, L, L]`` ``[…, 3, 3, 3]`` ============ ==================== ========================== For S3/S4 the projection is the tensor product of independent 1-D Fourier projections on each orientation axis:: result = scat_op.synthesis(xnorm, fft_ang=True, NUM_EPOCHS=300) * **fft_nharm** (:py:class:`int`, *optional*) -- Number of Fourier harmonics to keep beyond the DC term when ``fft_ang=True``. Default is ``1``. The number of output coefficients per orientation axis is ``1 + 2*fft_nharm`` (with ``fft_imaginary=True``) or ``1 + fft_nharm`` (with ``fft_imaginary=False``). * **fft_imaginary** (:py:class:`bool`, *optional*) -- Whether to include both cosine **and** sine components when ``fft_ang=True``. Default is ``True``. - ``True`` *(recommended)*: output per axis = ``1 + 2*fft_nharm`` coefficients ``[DC, cos₁, sin₁, cos₂, sin₂, …]``. The harmonic amplitude :math:`\sqrt{c_k^2 + s_k^2}` is independent of the image orientation. - ``False``: output per axis = ``1 + fft_nharm`` coefficients ``[DC, cos₁, cos₂, …]``. A field oriented at 90° (zero of cosine) gives a near-zero first-harmonic coefficient even if strongly anisotropic — use only if orientation is fixed. * **EVAL_FREQUENCY** (:py:class:`int`, *optional*) -- Print the current loss every N L-BFGS-B iterations. Default is ``100``. * **NUM_EPOCHS** (:py:class:`int`, *optional*) -- Maximum number of L-BFGS-B iterations **per resolution level**. Default is ``300``. The optimiser may stop earlier when its convergence criterion is met. Total wall-clock time scales as ``nstep × NUM_EPOCHS``. * **scat_cov_method** (:py:class:`str`, *optional*) -- Internal method for computing scattering covariances. Default is ``'eval'`` (recommended), which uses :meth:`eval` with ``norm='auto'`` and caches the normalisation after the first call. Any other value falls back to the legacy ``scattering_cov()`` path (2D only, kept for backward compatibility). * **n_up** (:py:class:`int`, *optional*) -- Number of **extra upsampling steps** beyond the target size, keeping the same ``Jmax``. Default is ``0`` (no extra steps). With ``n_up=1``, after completing synthesis at N × N, the algorithm continues at 2N × 2N using the same scattering statistics as the N × N target. The result is a field that locally matches the statistics of the original target, embedded in a larger canvas. Only available for 2D maps. +--------+-----------------------------+ | n_up | Output size (N × N target) | +========+=============================+ | 0 | N × N (standard) | +--------+-----------------------------+ | 1 | 2N × 2N | +--------+-----------------------------+ | 2 | 4N × 4N | +--------+-----------------------------+ The ``Jmax`` used during n_up steps is pinned to the value effective for the N × N target so that wavelet filters and the norm cache remain consistent:: result = scat_op.synthesis(xnorm_256, nstep=3, n_up=1) # result.shape == (512, 512) :returns: :py:class:`numpy.ndarray` -- Synthesised field. - Shape equals ``image_target.shape`` when ``synthesised_N=1`` and ``n_up=0``. - When ``synthesised_N > 1``, a leading dimension of size ``synthesised_N`` is added. - When ``n_up > 0`` (2D only), spatial dimensions are multiplied by ``2**n_up``. .. admonition:: Examples **Minimal 2D synthesis**:: import foscat.scat_cov2D as sc import numpy as np scat_op = sc.funct(NORIENT=4) xnorm = (image - image.mean()) / image.std() result = scat_op.synthesis(xnorm, seed=10, nstep=3, NUM_EPOCHS=300) **Batch synthesis with mask**:: mask = np.ones_like(xnorm) mask[invalid] = 0.0 results = scat_op.synthesis( xnorm, in_mask=mask, synthesised_N=4, nstep=3, iso_ang=True, NUM_EPOCHS=500, ) # results.shape == (4, H, W) **Inpainting — reconstruct missing pixels, freeze observed ones**:: grd_mask = np.zeros_like(xnorm) grd_mask[hole_pixels] = 1.0 result = scat_op.synthesis( xnorm, grd_mask=grd_mask, nstep=3, edge=True, NUM_EPOCHS=500 ) **Upsampled synthesis (n_up)**:: result = scat_op.synthesis(xnorm_256, nstep=3, n_up=1, NUM_EPOCHS=300) # result.shape == (512, 512) **Cross-covariance — component separation**:: result = scat_op.synthesis( cmb_estimate, reference=dust_template, nstep=3, iso_ang=True, NUM_EPOCHS=300, ) **Isotropic Gaussianised synthesis**:: result = scat_op.synthesis( image_target, to_gaussian=True, iso_ang=True, nstep=4, NUM_EPOCHS=500, ) .. admonition:: Notes The optimiser is ``scipy.optimize.fmin_l_bfgs_b``, a quasi-Newton method that uses a low-memory approximation of the inverse Hessian. It converges in tens to a few hundred iterations for typical scattering-covariance losses, far fewer than first-order methods. For large maps (nside ≥ 256 for HEALPix, H/W ≥ 256 for 2D) always use ``nstep ≥ 3``. Starting directly at full resolution wastes compute and converges poorly. .. py:method:: to_numpy(x)