13961110627

Committed 20 Mar 2025 03:06AM UTC coverage: 69.205% (-0.8%) from 70.039%

Build # 13961110627

Build Type

push

github

Committed by

avivajpeyi

Commit Message

increase the py version

Run Details

150 of 244 branches covered (61.48%)

Branch coverage included in aggregate %.

868 of 1227 relevant lines covered (70.74%)

0.71 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.74

/src/pywavelet/transforms/phi_computer.py

import numpy as np
from jaxtyping import Array, Float

from ..backend import betainc, ifft, xp

__all__ = ["phitilde_vec_norm", "phi_vec", "omega"]


def omega(Nf: int, Nt: int) -> Float[Array, " Nt//2+1"]:
    """Get the angular frequencies of the time domain signal."""
    df = 2 * np.pi / (Nf * Nt)
    return df * np.arange(0, Nt // 2 + 1)


def phitilde_vec_norm(Nf: int, Nt: int, d: float) -> Float[Array, " Nt//2+1"]:
    """Normalize phitilde for inverse frequency domain transform."""
    omegas = omega(Nf, Nt)
    _phi_t = _phitilde_vec(omegas, Nf, d) / (np.pi ** (-1 / 2))
    return xp.array(_phi_t)


def phi_vec(Nf: int, d: float = 4.0, q: int = 16) -> Float[Array, " 2*q*Nf"]:
    """get time domain phi as fourier transform of _phitilde_vec
    q: number of Nf bins over which the window extends?

    """
    insDOM = 1.0 / np.sqrt(np.pi / Nf)
    K = q * 2 * Nf
    half_K = q * Nf  # xp.int64(K/2)

    dom = 2 * np.pi / K  # max frequency is K/2*dom = pi/dt = OM

    DX = np.zeros(K, dtype=np.complex128)

    # zero frequency
    DX[0] = insDOM

    DX = DX.copy()
    # postive frequencies
    DX[1 : half_K + 1] = _phitilde_vec(dom * np.arange(1, half_K + 1), Nf, d)
    # negative frequencies
    DX[half_K + 1 :] = _phitilde_vec(
        -dom * np.arange(half_K - 1, 0, -1), Nf, d
    )
    DX = K * ifft(DX, K)

    phi = np.zeros(K)
    phi[0:half_K] = np.real(DX[half_K:K])
    phi[half_K:] = np.real(DX[0:half_K])

    nrm = np.sqrt(2.0) / np.sqrt(K / dom)  # *xp.linalg.norm(phi)

    phi *= nrm
    return xp.array(phi)


def _phitilde_vec(
    omega: Float[Array, " Nt//2+1"], Nf: int, d: float = 4.0
) -> Float[Array, " Nt//2+1"]:
    """Compute phi_tilde(omega_i) array, nx is filter steepness, defaults to 4.

    Eq 11 of https://arxiv.org/pdf/2009.00043.pdf (Cornish et al. 2020)

    phi(omega_i) =
        1/sqrt(2π∆F) if |omega_i| < A
        1/sqrt(2π∆F) cos(nu_d π/2 * |omega|-A / B) if A < |omega_i| < A + B

    Where nu_d = normalized incomplete beta function



    Parameters
    ----------
    ω : xp.ndarray
        Array of angular frequencies
    Nf : int
        Number of frequency bins
    d : float, optional
        Number of standard deviations for the gaussian wavelet, by default 4.

    Returns
    -------
    xp.ndarray
        Array of phi_tilde(omega_i) values

    """
    dF = 1.0 / (2 * Nf)  # NOTE: missing 1/dt?
    dOmega = 2 * np.pi * dF  # Near Eq 10 # 2 pi times DF
    inverse_sqrt_dOmega = 1.0 / np.sqrt(dOmega)

    A = dOmega / 4
    B = dOmega - 2 * A  # Cannot have B \leq 0.
    if B <= 0:
        raise ValueError("B must be greater than 0")

    phi = np.zeros(omega.size)
    mask = (A <= np.abs(omega)) & (np.abs(omega) < A + B)  # Minor changes
    vd = (np.pi / 2.0) * _nu_d(omega[mask], A, B, d=d)  # different from paper
    phi[mask] = inverse_sqrt_dOmega * xp.cos(vd)
    phi[np.abs(omega) < A] = inverse_sqrt_dOmega
    return phi


def _nu_d(
    omega: Float[Array, " Nt//2+1"], A: float, B: float, d: float = 4.0
) -> Float[Array, " Nt//2+1"]:
    """Compute the normalized incomplete beta function.

    Parameters
    ----------
    ω : xp.ndarray
        Array of angular frequencies
    A : float
        Lower bound for the beta function
    B : float
        Upper bound for the beta function
    d : float, optional
        Number of standard deviations for the gaussian wavelet, by default 4.

    Returns
    -------
    xp.ndarray
        Array of ν_d values

    scipy.special.betainc
    https://docs.scipy.org/doc/scipy-1.7.1/reference/reference/generated/scipy.special.betainc.html

    """
    x = (np.abs(omega) - A) / B
    return betainc(d, d, x)

1	import numpy as np	1✔
2	from jaxtyping import Array, Float	1✔
3
4	from ..backend import betainc, ifft, xp	1✔
5
6	__all__ = ["phitilde_vec_norm", "phi_vec", "omega"]	1✔
7
8
9	def omega(Nf: int, Nt: int) -> Float[Array, " Nt//2+1"]:	1✔
10	"""Get the angular frequencies of the time domain signal."""
11	df = 2 * np.pi / (Nf * Nt)	1✔
12	return df * np.arange(0, Nt // 2 + 1)	1✔
13
14
15	def phitilde_vec_norm(Nf: int, Nt: int, d: float) -> Float[Array, " Nt//2+1"]:	1✔
16	"""Normalize phitilde for inverse frequency domain transform."""
17	omegas = omega(Nf, Nt)	1✔
18	_phi_t = _phitilde_vec(omegas, Nf, d) / (np.pi ** (-1 / 2))	1✔
19	return xp.array(_phi_t)	1✔
20
21
22	def phi_vec(Nf: int, d: float = 4.0, q: int = 16) -> Float[Array, " 2qNf"]:	1✔
23	"""get time domain phi as fourier transform of _phitilde_vec
24	q: number of Nf bins over which the window extends?
25
26	"""
27	insDOM = 1.0 / np.sqrt(np.pi / Nf)	1✔
28	K = q * 2 * Nf	1✔
29	half_K = q * Nf # xp.int64(K/2)	1✔
30
31	dom = 2 * np.pi / K # max frequency is K/2*dom = pi/dt = OM	1✔
32
33	DX = np.zeros(K, dtype=np.complex128)	1✔
34
35	# zero frequency
36	DX[0] = insDOM	1✔
37
38	DX = DX.copy()	1✔
39	# postive frequencies
40	DX[1 : half_K + 1] = _phitilde_vec(dom * np.arange(1, half_K + 1), Nf, d)	1✔
41	# negative frequencies
42	DX[half_K + 1 :] = _phitilde_vec(	1✔
43	-dom * np.arange(half_K - 1, 0, -1), Nf, d
44	)
45	DX = K * ifft(DX, K)	1✔
46
47	phi = np.zeros(K)	1✔
48	phi[0:half_K] = np.real(DX[half_K:K])	1✔
49	phi[half_K:] = np.real(DX[0:half_K])	1✔
50
51	nrm = np.sqrt(2.0) / np.sqrt(K / dom) # *xp.linalg.norm(phi)	1✔
52
53	phi *= nrm	1✔
54	return xp.array(phi)	1✔
55
56
57	def _phitilde_vec(	1✔
58	omega: Float[Array, " Nt//2+1"], Nf: int, d: float = 4.0
59	) -> Float[Array, " Nt//2+1"]:
60	"""Compute phi_tilde(omega_i) array, nx is filter steepness, defaults to 4.
61
62	Eq 11 of https://arxiv.org/pdf/2009.00043.pdf (Cornish et al. 2020)
63
64	phi(omega_i) =
65	1/sqrt(2π∆F) if \|omega_i\| < A
66	1/sqrt(2π∆F) cos(nu_d π/2 * \|omega\|-A / B) if A < \|omega_i\| < A + B
67
68	Where nu_d = normalized incomplete beta function
69
70
71
72	Parameters
73	----------
74	ω : xp.ndarray
75	Array of angular frequencies
76	Nf : int
77	Number of frequency bins
78	d : float, optional
79	Number of standard deviations for the gaussian wavelet, by default 4.
80
81	Returns
82	-------
83	xp.ndarray
84	Array of phi_tilde(omega_i) values
85
86	"""
87	dF = 1.0 / (2 * Nf) # NOTE: missing 1/dt?	1✔
88	dOmega = 2 * np.pi * dF # Near Eq 10 # 2 pi times DF	1✔
89	inverse_sqrt_dOmega = 1.0 / np.sqrt(dOmega)	1✔
90
91	A = dOmega / 4	1✔
92	B = dOmega - 2 * A # Cannot have B \leq 0.	1✔
93	if B <= 0:	1!
94	raise ValueError("B must be greater than 0")	×
95
96	phi = np.zeros(omega.size)	1✔
97	mask = (A <= np.abs(omega)) & (np.abs(omega) < A + B) # Minor changes	1✔
98	vd = (np.pi / 2.0) * _nu_d(omega[mask], A, B, d=d) # different from paper	1✔
99	phi[mask] = inverse_sqrt_dOmega * xp.cos(vd)	1✔
100	phi[np.abs(omega) < A] = inverse_sqrt_dOmega	1✔
101	return phi	1✔
102
103
104	def _nu_d(	1✔
105	omega: Float[Array, " Nt//2+1"], A: float, B: float, d: float = 4.0
106	) -> Float[Array, " Nt//2+1"]:
107	"""Compute the normalized incomplete beta function.
108
109	Parameters
110	----------
111	ω : xp.ndarray
112	Array of angular frequencies
113	A : float
114	Lower bound for the beta function
115	B : float
116	Upper bound for the beta function
117	d : float, optional
118	Number of standard deviations for the gaussian wavelet, by default 4.
119
120	Returns
121	-------
122	xp.ndarray
123	Array of ν_d values
124
125	scipy.special.betainc
126	https://docs.scipy.org/doc/scipy-1.7.1/reference/reference/generated/scipy.special.betainc.html
127
128	"""
129	x = (np.abs(omega) - A) / B	1✔
130	return betainc(d, d, x)	1✔

avivajpeyi / pywavelet / 13961110627

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous