Utility Functions

VOICEBOX configuration, filesystem helpers, numeric routines, and combinatorics.

VOICEBOX configuration and environment

v_voicebox

V_VOICEBOX - Global parameters for Voicebox functions.

Set/get global parameters used by other functions in the VOICEBOX toolbox.

v_voicebox

v_voicebox(field=None, value=None) -> Any

Get or set global Voicebox parameters.

Parameters:

Name	Type	Description	Default
field	str	Parameter name. If None, returns all parameters.	None
value	optional	New value for the parameter.	None

Returns:

Type	Description
Result depends on inputs:	No args: returns dict of all parameters field only: returns value of that field (None if not found) field and value: sets field and returns all parameters

Source code in pyvoicebox/v_voicebox.py

def v_voicebox(field=None, value=None) -> Any:
    """Get or set global Voicebox parameters.

    Parameters
    ----------
    field : str, optional
        Parameter name. If None, returns all parameters.
    value : optional
        New value for the parameter.

    Returns
    -------
    Result depends on inputs:
        - No args: returns dict of all parameters
        - field only: returns value of that field (None if not found)
        - field and value: sets field and returns all parameters
    """
    global _PP
    if _PP is None:
        _PP = _init_defaults()

    if field is None:
        return dict(_PP)
    elif value is None:
        return _PP.get(field, None)
    else:
        if field not in _PP:
            raise ValueError(f"'{field}' is not a valid voicebox field name")
        _PP[field] = value
        return dict(_PP)

v_voicebox_update

V_VOICEBOX_UPDATE - Check for voicebox updates (stub).

v_voicebox_update

v_voicebox_update() -> None

Check for voicebox updates.

This is a stub. The MATLAB version checks for toolbox updates. In the Python port, use pip for package management.

Raises:

Type	Description
NotImplementedError	Always raised as this is a MATLAB-specific function.

Source code in pyvoicebox/v_voicebox_update.py

def v_voicebox_update() -> None:
    """Check for voicebox updates.

    This is a stub. The MATLAB version checks for toolbox updates.
    In the Python port, use pip for package management.

    Raises
    ------
    NotImplementedError
        Always raised as this is a MATLAB-specific function.
    """
    raise NotImplementedError(
        "v_voicebox_update is MATLAB-specific. Use pip to manage Python packages."
    )

v_paramsetch

V_PARAMSETCH - Set parameters for speech processing algorithms (stub).

v_paramsetch

v_paramsetch(*args, **kwargs) -> dict

Set parameters for speech processing algorithms.

The full MATLAB implementation handles parameter sets for various speech processing configurations. A stub is provided.

Raises:

Type	Description
NotImplementedError	Full implementation pending.

Source code in pyvoicebox/v_paramsetch.py

def v_paramsetch(*args, **kwargs) -> dict:
    """Set parameters for speech processing algorithms.

    The full MATLAB implementation handles parameter sets for various
    speech processing configurations. A stub is provided.

    Raises
    ------
    NotImplementedError
        Full implementation pending.
    """
    raise NotImplementedError(
        "v_paramsetch full implementation pending."
    )

v_hostipinfo

V_HOSTIPINFO - Get host name and IP info using Python equivalents.

v_hostipinfo

v_hostipinfo(m='') -> dict

Get host name and IP address info.

Parameters:

Name	Type	Description	Default
m	str	Mode string. 'h' for hostname, 'i' for IP address.	''

Returns:

Name	Type	Description
info	str	Hostname or IP address.

Source code in pyvoicebox/v_hostipinfo.py

def v_hostipinfo(m='') -> dict:
    """Get host name and IP address info.

    Parameters
    ----------
    m : str, optional
        Mode string. 'h' for hostname, 'i' for IP address.

    Returns
    -------
    info : str
        Hostname or IP address.
    """
    if not m or 'h' in m:
        return socket.gethostname()
    if 'i' in m:
        try:
            return socket.gethostbyname(socket.gethostname())
        except socket.gaierror:
            return '127.0.0.1'
    return socket.gethostname()

v_winenvar

V_WINENVAR - Read Windows environment variable (stub).

v_winenvar

v_winenvar(name) -> str

Read a Windows environment variable.

This is a cross-platform stub that uses os.environ.

Parameters:

Name	Type	Description	Default
name	str	Environment variable name.	required

Returns:

Name	Type	Description
value	str or None	Value of the environment variable, or None if not found.

Source code in pyvoicebox/v_winenvar.py

def v_winenvar(name) -> str:
    """Read a Windows environment variable.

    This is a cross-platform stub that uses os.environ.

    Parameters
    ----------
    name : str
        Environment variable name.

    Returns
    -------
    value : str or None
        Value of the environment variable, or None if not found.
    """
    return os.environ.get(name, None)

v_unixwhich

V_UNIXWHICH - Search system path for an executable (Python equivalent).

v_unixwhich

v_unixwhich(c, e=None) -> str

Search system path for an executable program.

Uses Python's shutil.which() as a cross-platform equivalent.

Parameters:

Name	Type	Description	Default
c	str	Name of file to search for.	required
e	str	Extensions to search (ignored; shutil.which handles this).	None

Returns:

Name	Type	Description
f	str or None	Full pathname of executable, or None if not found.

Source code in pyvoicebox/v_unixwhich.py

def v_unixwhich(c, e=None) -> str:
    """Search system path for an executable program.

    Uses Python's shutil.which() as a cross-platform equivalent.

    Parameters
    ----------
    c : str
        Name of file to search for.
    e : str, optional
        Extensions to search (ignored; shutil.which handles this).

    Returns
    -------
    f : str or None
        Full pathname of executable, or None if not found.
    """
    return shutil.which(c)

v_regexfiles

V_REGEXFILES - Find files matching a regular expression pattern.

v_regexfiles

v_regexfiles(
    pattern, directory=".", recursive=False
) -> list

Find files matching a regular expression pattern.

Parameters:

Name	Type	Description	Default
pattern	str	Regular expression pattern to match filenames.	required
directory	str	Directory to search. Default is current directory.	'.'
recursive	bool	If True, search recursively. Default False.	False

Returns:

Name	Type	Description
files	list of str	List of matching file paths.

Source code in pyvoicebox/v_regexfiles.py

def v_regexfiles(pattern, directory='.', recursive=False) -> list:
    """Find files matching a regular expression pattern.

    Parameters
    ----------
    pattern : str
        Regular expression pattern to match filenames.
    directory : str, optional
        Directory to search. Default is current directory.
    recursive : bool, optional
        If True, search recursively. Default False.

    Returns
    -------
    files : list of str
        List of matching file paths.
    """
    regex = re.compile(pattern)
    matches = []

    if recursive:
        for root, dirs, filenames in os.walk(directory):
            for fn in filenames:
                if regex.search(fn):
                    matches.append(os.path.join(root, fn))
    else:
        for fn in os.listdir(directory):
            if regex.search(fn) and os.path.isfile(os.path.join(directory, fn)):
                matches.append(os.path.join(directory, fn))

    return sorted(matches)

v_fopenmkd

V_FOPENMKD - Open file, creating directories if needed.

v_fopenmkd

v_fopenmkd(filename, mode='r', **kwargs) -> Any

Open a file, creating any missing parent directories.

Parameters:

Name	Type	Description	Default
filename	str	Path to file.	required
mode	str	File open mode (default 'r').	'r'
**kwargs		Additional arguments passed to open().	{}

Returns:

Name	Type	Description
fid	file object	Open file handle.

Source code in pyvoicebox/v_fopenmkd.py

def v_fopenmkd(filename, mode='r', **kwargs) -> Any:
    """Open a file, creating any missing parent directories.

    Parameters
    ----------
    filename : str
        Path to file.
    mode : str, optional
        File open mode (default 'r').
    **kwargs
        Additional arguments passed to open().

    Returns
    -------
    fid : file object
        Open file handle.
    """
    try:
        return open(filename, mode, **kwargs)
    except FileNotFoundError:
        parent = os.path.dirname(filename)
        if parent:
            os.makedirs(parent, exist_ok=True)
        return open(filename, mode, **kwargs)

v_finishat

V_FINISHAT - Print estimated finish time of a long computation.

v_finishat

v_finishat(frac, tol=None, fmt=None) -> str

Print estimated finish time of a long computation.

Parameters:

Name	Type	Description	Default
frac	float or ndarray	Fraction of total computation completed (0 to 1). 0 initializes the timer.	required
tol	float	Tolerance in minutes before reprinting.	None
fmt	str	Format string.	None

Returns:

Name	Type	Description
eta	str	Estimated finish time string.

Source code in pyvoicebox/v_finishat.py

def v_finishat(frac, tol=None, fmt=None) -> str:
    """Print estimated finish time of a long computation.

    Parameters
    ----------
    frac : float or ndarray
        Fraction of total computation completed (0 to 1).
        0 initializes the timer.
    tol : float, optional
        Tolerance in minutes before reprinting.
    fmt : str, optional
        Format string.

    Returns
    -------
    eta : str
        Estimated finish time string.
    """
    global _state

    if fmt is None:
        fmt = 'Estimated finish at {time} ({pct:.0f}% done, {remaining} remaining)'

    if frac <= 0:
        _state = {
            'tstart': time.time(),
            'oldt': 0,
            'oldnw': 0,
        }
        return 'Unknown'

    if 'tstart' not in _state:
        _state = {
            'tstart': time.time(),
            'oldt': 0,
            'oldnw': 0,
        }

    elapsed = time.time() - _state['tstart']
    if frac > 0:
        sectogo = (1.0 / frac - 1) * elapsed
    else:
        return 'Unknown'

    now = time.time()
    finish_time = datetime.now() + timedelta(seconds=sectogo)

    if tol is None:
        tol = max(0.1 * sectogo / 60, 1)

    oldt = _state.get('oldt', 0)
    oldnw = _state.get('oldnw', 0)

    if (oldt == 0 or
            abs(finish_time.timestamp() - oldt) > tol * 60 and (now - oldnw) > 10 or
            (now - oldnw) > 600):

        _state['oldt'] = finish_time.timestamp()
        _state['oldnw'] = now

        if finish_time.date() == datetime.now().date():
            eta = finish_time.strftime('%H:%M')
        else:
            eta = finish_time.strftime('%H:%M %d-%b-%Y')

        # Format remaining time
        if sectogo >= 3600:
            remaining = f'{sectogo / 3600:.1f} hr'
        elif sectogo >= 60:
            remaining = f'{sectogo / 60:.1f} min'
        else:
            remaining = f'{sectogo:.0f} sec'

        msg = fmt.format(time=eta, pct=frac * 100, remaining=remaining)
        print(msg, file=sys.stderr)
        return eta

    return ''

v_m2htmlpwd

V_M2HTMLPWD - MATLAB-specific HTML documentation utility (stub).

v_m2htmlpwd

v_m2htmlpwd() -> str

Generate HTML documentation from MATLAB m-files.

This is a MATLAB-specific function and has no Python equivalent.

Raises:

Type	Description
NotImplementedError	Always raised as this is MATLAB-specific.

Source code in pyvoicebox/v_m2htmlpwd.py

def v_m2htmlpwd() -> str:
    """Generate HTML documentation from MATLAB m-files.

    This is a MATLAB-specific function and has no Python equivalent.

    Raises
    ------
    NotImplementedError
        Always raised as this is MATLAB-specific.
    """
    raise NotImplementedError(
        "v_m2htmlpwd is MATLAB-specific. Use Sphinx or pydoc for Python documentation."
    )

Numerical helpers

v_atan2sc

V_ATAN2SC - sin and cosine of atan(y/x).

v_atan2sc

v_atan2sc(
    y, x
) -> tuple[ndarray, ndarray, ndarray, ndarray]

Compute sin and cosine of atan(y/x) robustly.

Parameters:

Name	Type	Description	Default
y	array_like	Input values (same shape).	required
x	array_like	Input values (same shape).	required

Returns:

Name	Type	Description
s	ndarray	sin(t) where tan(t) = y/x
c	ndarray	cos(t) where tan(t) = y/x
r	ndarray	sqrt(x^2 + y^2)
t	ndarray	atan2(y, x)

Source code in pyvoicebox/v_atan2sc.py

def v_atan2sc(y, x) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Compute sin and cosine of atan(y/x) robustly.

    Parameters
    ----------
    y, x : array_like
        Input values (same shape).

    Returns
    -------
    s : ndarray
        sin(t) where tan(t) = y/x
    c : ndarray
        cos(t) where tan(t) = y/x
    r : ndarray
        sqrt(x^2 + y^2)
    t : ndarray
        atan2(y, x)
    """
    y = np.asarray(y, dtype=float)
    x = np.asarray(x, dtype=float)

    s = np.zeros_like(y)
    c = np.full_like(y, np.nan)
    r = np.zeros_like(y)
    t = np.full_like(y, np.nan)

    # Handle y == 0
    m = (y == 0)
    if np.any(m):
        neg_x = (x[m] < 0).astype(float)
        t[m] = neg_x * np.pi
        c[m] = 1.0 - 2.0 * neg_x
        r[m] = np.abs(x[m])

    # Handle |y| > |x| and not yet computed
    m2 = (np.abs(y) > np.abs(x)) & np.isnan(c)
    if np.any(m2):
        q = x[m2] / y[m2]
        u = np.sqrt(1.0 + q ** 2) * np.sign(y[m2])
        s[m2] = 1.0 / u
        c[m2] = s[m2] * q
        r[m2] = y[m2] * u

    # Handle remaining (|x| >= |y|, y != 0)
    m3 = np.isnan(c)
    if np.any(m3):
        q = y[m3] / x[m3]
        u = np.sqrt(1.0 + q ** 2) * np.sign(x[m3])
        c[m3] = 1.0 / u
        s[m3] = c[m3] * q
        r[m3] = x[m3] * u

    # Compute t where still NaN
    m4 = np.isnan(t)
    if np.any(m4):
        t[m4] = np.arctan2(s[m4], c[m4])

    return s, c, r, t

v_logsum

V_LOGSUM - log(sum(k.*exp(x),d)) computed avoiding overflow/underflow.

v_logsum

v_logsum(x, d=None, k=None) -> ndarray

Compute log(sum(k.*exp(x), d)) avoiding overflow/underflow.

Parameters:

Name	Type	Description	Default
x	array_like	Data array.	required
d	int	Axis to sum along. Default: first non-singleton dimension.	None
k	array_like	Scaling array (same shape as x, or vector along axis d).	None

Returns:

Name	Type	Description
y	ndarray	log(sum(k.*exp(x), d))

Source code in pyvoicebox/v_logsum.py

def v_logsum(x, d=None, k=None) -> np.ndarray:
    """Compute log(sum(k.*exp(x), d)) avoiding overflow/underflow.

    Parameters
    ----------
    x : array_like
        Data array.
    d : int, optional
        Axis to sum along. Default: first non-singleton dimension.
    k : array_like, optional
        Scaling array (same shape as x, or vector along axis d).

    Returns
    -------
    y : ndarray
        log(sum(k.*exp(x), d))
    """
    x = np.asarray(x, dtype=float)

    if d is None:
        d = first_nonsingleton(x)

    n = x.shape[d]
    if n <= 1:
        if k is None:
            return x.copy()
        else:
            return x + np.log(np.asarray(k, dtype=float))

    q = np.max(x, axis=d, keepdims=True)

    # Handle -inf slices: replace -inf max with 0 for safe subtraction,
    # then restore -inf in the result. This avoids (-inf) - (-inf) = nan.
    a = np.isinf(q)
    q_safe = np.where(a, 0.0, q)

    if k is None:
        y = q_safe + np.log(np.sum(np.exp(x - q_safe), axis=d, keepdims=True))
    else:
        k = np.asarray(k, dtype=float)
        # Broadcast k to match x along axis d
        if k.ndim == 1 and k.shape[0] == n:
            shape = [1] * x.ndim
            shape[d] = n
            k = k.reshape(shape)
        y = q_safe + np.log(np.sum(np.exp(x - q_safe) * k, axis=d, keepdims=True))

    # Restore +-inf for slices that had infinite max
    y = np.where(a, q, y)

    # Remove the keepdims axis to match MATLAB behavior
    y = np.squeeze(y, axis=d)

    return y

v_gammalns

V_GAMMALNS - Log of Gamma(x) for positive or negative real x.

v_gammalns

v_gammalns(x, return_sign=False) -> ndarray

Compute log(|Gamma(x)|) and optionally sign(Gamma(x)).

Parameters:

Name	Type	Description	Default
x	array_like	Real input values.	required
return_sign	bool	If True, return (y, s) where s = sign(Gamma(x)). If False (default), return y which may be complex for negative Gamma.	False

Returns:

Name	Type	Description
y	ndarray	log(|Gamma(x)|) if return_sign=True, else log(Gamma(x)) (complex where needed).
s	ndarray (only if return_sign=True)	sign(Gamma(x))

Source code in pyvoicebox/v_gammalns.py

def v_gammalns(x, return_sign=False) -> np.ndarray:
    """Compute log(|Gamma(x)|) and optionally sign(Gamma(x)).

    Parameters
    ----------
    x : array_like
        Real input values.
    return_sign : bool, optional
        If True, return (y, s) where s = sign(Gamma(x)).
        If False (default), return y which may be complex for negative Gamma.

    Returns
    -------
    y : ndarray
        log(|Gamma(x)|) if return_sign=True, else log(Gamma(x)) (complex where needed).
    s : ndarray (only if return_sign=True)
        sign(Gamma(x))
    """
    x = np.asarray(x, dtype=float)
    y = np.zeros_like(x)
    s = np.ones_like(x)

    m = x <= 0
    # Non-negative x values
    if not np.all(m):
        y[~m] = gammaln(x[~m])

    # Non-positive integers: Gamma is infinite
    f = m & (x == np.fix(x))
    if np.any(f):
        y[f] = np.inf
        m = m & ~f

    # Negative non-integer x values
    if np.any(m):
        t = np.sin(np.pi * x[m])
        if return_sign:
            p = t < 0
            s[m] = 1 - 2 * p
            y[m] = np.log(np.pi) - gammaln(1 - x[m]) - np.log(np.abs(t))
        else:
            y[m] = np.log(np.pi) - gammaln(1 - x[m]) - np.log(t + 0j)
            # Make result complex where needed
            y = y.astype(complex)

    if return_sign:
        return y, s
    return y

v_hypergeom1f1

V_HYPERGEOM1F1 - Confluent hypergeometric function 1F1 (Kummer's M).

v_hypergeom1f1

v_hypergeom1f1(
    a, b, z, tol=1e-10, maxj=500, th=30
) -> tuple[ndarray, ndarray]

Confluent hypergeometric function 1F1(a; b; z) = M(a; b; z).

Parameters:

Name	Type	Description	Default
a	float	First parameter (real scalar).	required
b	float	Second parameter (real scalar).	required
z	array_like	Input values (real).	required
tol	float	Tolerance. Default 1e-10.	1e-10
maxj	int	Maximum iterations. Default 500.	500
th	float	Threshold for algorithm selection. Default 30.	30

Returns:

Name	Type	Description
h	ndarray	Output values.
l	ndarray	Number of iterations used.

Source code in pyvoicebox/v_hypergeom1f1.py

def v_hypergeom1f1(a, b, z, tol=1e-10, maxj=500, th=30) -> tuple[np.ndarray, np.ndarray]:
    """Confluent hypergeometric function 1F1(a; b; z) = M(a; b; z).

    Parameters
    ----------
    a : float
        First parameter (real scalar).
    b : float
        Second parameter (real scalar).
    z : array_like
        Input values (real).
    tol : float, optional
        Tolerance. Default 1e-10.
    maxj : int, optional
        Maximum iterations. Default 500.
    th : float, optional
        Threshold for algorithm selection. Default 30.

    Returns
    -------
    h : ndarray
        Output values.
    l : ndarray
        Number of iterations used.
    """
    z = np.asarray(z, dtype=float)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)
    h = np.zeros_like(z, dtype=float)
    l = np.zeros_like(z, dtype=int)

    a1 = a - 1
    b1 = b - 1
    ba = b - a

    for idx in np.ndindex(z.shape):
        y = z[idx]
        q = False  # break criterion

        if abs(y) < th:
            # Algorithm 1: series expansion (13.2.2)
            d = 1.0
            g = 1.0
            jlim = 0
            k = (b1 + y) ** 2 - 4 * a1 * y
            if k >= 0:
                jlim = max(jlim, 0.5 * (np.sqrt(k) - (b1 + y)))
            k = (b1 - y) ** 2 + 4 * a1 * y
            if k >= 0:
                jlim = max(jlim, 0.5 * (np.sqrt(k) - (b1 - y)))
            jlim = min(maxj, jlim)

            for j in range(1, maxj + 1):
                d = d * y * (a1 + j) / (j * (b1 + j))
                g = g + d
                p = abs(d) < tol * abs(g)
                if q and p and j >= jlim:
                    break
                q = p

        elif y > 0:
            # Algorithm 2: asymptotic for large positive y
            d = 1.0
            g = 1.0
            jlim = 1
            k = (ba - 1 - a + y) ** 2 + 4 * a * (ba - 1)
            if k >= 0:
                jlim = max(jlim, 0.5 * (np.sqrt(k) - (ba - a - 1 + y)))
            k = (ba - a - 1 - y) ** 2 + 4 * a * (ba - 1)
            if k >= 0:
                jlim = max(jlim, 0.5 * (np.sqrt(k) - (ba - a - 1 - y)))
            jlim = int(min(maxj, jlim))

            for j in range(1, jlim + 1):
                d = d * (ba - 1 + j) * (j - a) / (j * y)
                g = g + d
                p = abs(d) < tol * abs(g)
                if q and p:
                    break
                q = p

            gl, gs = v_gammalns(np.array([a, b]))
            g = gs[0] * gs[1] * np.exp(y + gl[1] - gl[0] + (a - b) * np.log(y)) * g

        else:
            # Algorithm 3: large negative y using M(a;b;z) = exp(z)*M(b-a;b;-z)
            d = 1.0
            g = 1.0
            jlim = 1
            k = (a1 - ba + y) ** 2 + 4 * a1 * ba
            if k >= 0:
                jlim = max(jlim, 0.5 * (np.sqrt(k) - (a1 - ba + y)))
            k = (a1 - ba - y) ** 2 + 4 * a1 * ba
            if k >= 0:
                jlim = max(jlim, 0.5 * (np.sqrt(k) - (a1 - ba - y)))
            jlim = int(min(maxj, jlim))

            for j in range(1, jlim + 1):
                d = d * (a1 + j) * (ba - j) / (j * y)
                g = g + d
                p = abs(d) < tol * abs(g)
                if q and p:
                    break
                q = p

            gl, gs = v_gammalns(np.array([ba, b]))
            g = gs[0] * gs[1] * np.exp(gl[1] - gl[0] - a * np.log(-y)) * g

        h[idx] = g
        l[idx] = j

    if scalar_input:
        return h.item(), l.item()
    return h, l

v_dualdiag

V_DUALDIAG - Simultaneous diagonalization of two Hermitian matrices.

v_dualdiag

v_dualdiag(w, b) -> tuple[ndarray, ndarray, ndarray]

Simultaneous diagonalization of two Hermitian matrices.

Parameters:

Name	Type	Description	Default
w	(array_like, shape(n, n))	Hermitian matrix.	required
b	(array_like, shape(n, n))	Hermitian matrix.	required

Returns:

Name	Type	Description
a	(ndarray, shape(n, n))	Diagonalizing matrix.
d	(ndarray, shape(n))	Real diagonal elements: A'BA = diag(D).
e	(ndarray, shape(n))	Real diagonal elements: A'WA = diag(E).

Source code in pyvoicebox/v_dualdiag.py

def v_dualdiag(w, b) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Simultaneous diagonalization of two Hermitian matrices.

    Parameters
    ----------
    w : array_like, shape (n, n)
        Hermitian matrix.
    b : array_like, shape (n, n)
        Hermitian matrix.

    Returns
    -------
    a : ndarray, shape (n, n)
        Diagonalizing matrix.
    d : ndarray, shape (n,)
        Real diagonal elements: A'*B*A = diag(D).
    e : ndarray, shape (n,)
        Real diagonal elements: A'*W*A = diag(E).
    """
    w = np.asarray(w, dtype=complex if np.iscomplexobj(w) else float)
    b = np.asarray(b, dtype=complex if np.iscomplexobj(b) else float)
    n = w.shape[0]

    # Generalized eigendecomposition
    eigenvalues, eigenvectors = eig(b + b.conj().T, w + w.conj().T)

    a = eigenvectors

    if np.isrealobj(eigenvalues):
        eigenvalues = eigenvalues.real
        # Sort by absolute value descending
        idx = np.argsort(np.abs(eigenvalues))[::-1]
        a = a[:, idx]
        e = np.real(np.diag(a.conj().T @ w @ a))
        d = np.real(np.diag(a.conj().T @ b @ a))
    else:
        d = a.conj().T @ b @ a
        e = a.conj().T @ w @ a

    return a, d, e

v_mintrace

V_MINTRACE - Find row permutation to minimize trace.

v_mintrace

v_mintrace(x) -> ndarray

Find row permutation to minimize trace of x(p,:).

Parameters:

Name	Type	Description	Default
x	(array_like, shape(n, n))	Square real-valued matrix.	required

Returns:

Name	Type	Description
p	(ndarray, shape(n))	Row permutation that minimizes trace(x[p, :]).

Source code in pyvoicebox/v_mintrace.py

def v_mintrace(x) -> np.ndarray:
    """Find row permutation to minimize trace of x(p,:).

    Parameters
    ----------
    x : array_like, shape (n, n)
        Square real-valued matrix.

    Returns
    -------
    p : ndarray, shape (n,)
        Row permutation that minimizes trace(x[p, :]).
    """
    x = np.asarray(x, dtype=float)
    n = x.shape[0]
    perms = v_permutes(n)  # shape (n!, n)
    # Compute trace for each permutation
    cols = np.arange(n)
    traces = np.sum(x[perms, cols], axis=1)
    best = np.argmin(traces)
    return perms[best, :]

v_quadpeak

V_QUADPEAK - Find quadratically-interpolated peak in an N-D array.

v_quadpeak

v_quadpeak(z) -> tuple[float, ndarray, int, ndarray]

Find quadratically-interpolated peak in an N-D array.

Parameters:

Name	Type	Description	Default
z	array_like	Input array (at least 3 elements in each non-singleton dimension).	required

Returns:

Name	Type	Description
v	float	Peak value.
x	ndarray	Position of peak (fractional subscript values, 1-based).
t	int	-1 for maximum, 0 for saddle point, +1 for minimum.
m	ndarray	Matrix defining the fitted quadratic.

Source code in pyvoicebox/v_quadpeak.py

def v_quadpeak(z) -> tuple[float, np.ndarray, int, np.ndarray]:
    """Find quadratically-interpolated peak in an N-D array.

    Parameters
    ----------
    z : array_like
        Input array (at least 3 elements in each non-singleton dimension).

    Returns
    -------
    v : float
        Peak value.
    x : ndarray
        Position of peak (fractional subscript values, 1-based).
    t : int
        -1 for maximum, 0 for saddle point, +1 for minimum.
    m : ndarray
        Matrix defining the fitted quadratic.
    """
    z = np.asarray(z, dtype=float)
    sz_orig = z.shape
    dz_orig = len(sz_orig)

    # In MATLAB, a row vector [1 3 5 4 2] has size [1, 5].
    # In Python, it has shape (5,). Treat 1D arrays as row vectors.
    if dz_orig == 1:
        sz = (1, sz_orig[0])
    else:
        sz = sz_orig
    dz = len(sz)
    psz = z.size

    # Find non-singleton dimensions (0-based)
    mz = [i for i in range(dz) if sz[i] > 1]
    nm = len(mz)
    vz = [sz[i] for i in mz]
    dx = max(mz) + 1 if mz else 1

    if nm == 0:
        raise ValueError('Cannot find peak of a scalar')
    if min(vz) < 3:
        raise ValueError('Need at least 3 points in each non-singleton dimension')

    nc = (nm + 1) * (nm + 2) // 2

    # Build the A matrix using column-major (Fortran) index decomposition
    # to match MATLAB's linear indexing
    a_mat = np.ones((psz, nc))
    ix = np.arange(psz, dtype=float)

    for i in range(nm):
        i1 = i + 1  # MATLAB 1-based loop variable
        dim_size = sz[mz[i]]
        jx = np.floor(ix / dim_size)

        # Linear term column (0-based)
        lin_col = i + nc - nm - 1
        a_mat[:, lin_col] = 1 + ix - jx * dim_size
        ix = jx

        # Quadratic/cross term columns (0-based)
        col_start = (i1 * i1 - i1 + 2) // 2 - 1
        col_end = i1 * (i1 + 1) // 2 - 1  # inclusive

        # Linear terms from dimension 0..i (columns nc-nm-1 to lin_col inclusive)
        lin_range_start = nc - nm - 1
        lin_terms = a_mat[:, lin_range_start:lin_col + 1]  # shape (psz, i1)
        current_lin = a_mat[:, lin_col]  # shape (psz,)

        a_mat[:, col_start:col_end + 1] = lin_terms * current_lin[:, np.newaxis]

    # Solve for polynomial coefficients
    a_mat = np.linalg.solve(a_mat.T @ a_mat, a_mat.T)

    # Use Fortran-order ravel to match MATLAB's column-major linear indexing
    if dz_orig == 1:
        z_flat = z.ravel()
    else:
        z_flat = z.ravel(order='F')

    c = a_mat @ z_flat
    w = np.zeros((nm + 1, nm + 1))

    # Fill w matrix using MATLAB's column-major linear indexing
    idx = np.arange(1, nc + 1)
    j_vals = np.floor((np.sqrt(8 * idx - 7) - 1) / 2).astype(int)
    for k in range(nc):
        lin_idx = int(idx[k] + j_vals[k] * (2 * nm + 1 - j_vals[k]) / 2) - 1
        # MATLAB column-major indexing into (nm+1) x (nm+1) matrix
        row = lin_idx % (nm + 1)
        col = lin_idx // (nm + 1)
        w[row, col] = c[k]
    w = (w + w.T) / 2.0

    mr = w[:nm, :nm]
    we = w[:nm, nm]
    y = -np.linalg.solve(mr, we)
    v = float(y @ we + w[nm, nm])

    # Insert singleton dimensions
    x = np.zeros(dx)
    for i, mi in enumerate(mz):
        x[mi] = y[i]

    m_out = np.zeros((dx + 1, dx + 1))
    mz_ext = list(mz) + [dx]
    for i, mi in enumerate(mz_ext):
        for j_idx, mj in enumerate(mz_ext):
            m_out[mi, mj] = w[i, j_idx]

    ev = np.linalg.eigvalsh(mr)
    t = int(np.all(ev > 0)) - int(np.all(ev < 0))

    return v, x, t, m_out

v_peak2dquad

V_PEAK2DQUAD - Find quadratically-interpolated peak in a 2D array.

v_peak2dquad

v_peak2dquad(z) -> float

Find quadratically-interpolated peak in a 2D array.

This is a wrapper around v_quadpeak.

Parameters:

Name	Type	Description	Default
z	(array_like, shape(m, n))	Input 2D array.	required

Returns:

Name	Type	Description
v	float	Peak value.
xy	ndarray	Position of peak.
t	int	-1 for maximum, 0 for saddle point, +1 for minimum.
m	ndarray	Fitted quadratic matrix.

Source code in pyvoicebox/v_peak2dquad.py

def v_peak2dquad(z) -> float:
    """Find quadratically-interpolated peak in a 2D array.

    This is a wrapper around v_quadpeak.

    Parameters
    ----------
    z : array_like, shape (m, n)
        Input 2D array.

    Returns
    -------
    v : float
        Peak value.
    xy : ndarray
        Position of peak.
    t : int
        -1 for maximum, 0 for saddle point, +1 for minimum.
    m : ndarray
        Fitted quadratic matrix.
    """
    return v_quadpeak(z)

Combinatorics and sorting

v_choosenk

V_CHOOSENK - All choices of K elements from 0:N-1.

v_choosenk

v_choosenk(n, k) -> ndarray

Generate all choices of K elements from 0:N-1 in lexical order.

Note: Returns 0-based indices (unlike MATLAB's 1-based).

Parameters:

Name	Type	Description	Default
n	int	Range of elements (0 to n-1).	required
k	int	Number of elements to choose.	required

Returns:

Name	Type	Description
x	ndarray of shape (C(n,k), k)	Each row is a combination.

Source code in pyvoicebox/v_choosenk.py

def v_choosenk(n, k) -> np.ndarray:
    """Generate all choices of K elements from 0:N-1 in lexical order.

    Note: Returns 0-based indices (unlike MATLAB's 1-based).

    Parameters
    ----------
    n : int
        Range of elements (0 to n-1).
    k : int
        Number of elements to choose.

    Returns
    -------
    x : ndarray of shape (C(n,k), k)
        Each row is a combination.
    """
    if k > n:
        return np.array([]).reshape(0, max(k, 0))
    if k == 0:
        return np.array([]).reshape(1, 0)

    result = list(combinations(range(n), k))
    x = np.array(result, dtype=int)
    if x.ndim == 1:
        x = x.reshape(1, -1)
    return x

v_choosrnk

V_CHOOSRNK - All choices of K elements from 0:N-1 with replacement.

v_choosrnk

v_choosrnk(n, k) -> ndarray

Generate all choices of K elements from 0:N-1 with replacement.

Note: Returns 0-based indices (unlike MATLAB's 1-based).

Parameters:

Name	Type	Description	Default
n	int	Range of elements (0 to n-1).	required
k	int	Number of elements to choose.	required

Returns:

Name	Type	Description
x	ndarray	Each row is a combination with replacement.

Source code in pyvoicebox/v_choosrnk.py

def v_choosrnk(n, k) -> np.ndarray:
    """Generate all choices of K elements from 0:N-1 with replacement.

    Note: Returns 0-based indices (unlike MATLAB's 1-based).

    Parameters
    ----------
    n : int
        Range of elements (0 to n-1).
    k : int
        Number of elements to choose.

    Returns
    -------
    x : ndarray
        Each row is a combination with replacement.
    """
    x = v_choosenk(n + k - 1, k)
    x = x - np.arange(k)
    return x

v_permutes

V_PERMUTES - All N! permutations of 0:N-1 + signatures.

v_permutes

v_permutes(n, return_sign=False) -> ndarray

Generate all N! permutations of 0:N-1 in lexical order.

Note: Returns 0-based indices (unlike MATLAB's 1-based).

Parameters:

Name	Type	Description	Default
n	int	Number of elements to permute.	required
return_sign	bool	If True, also return the signature of each permutation.	False

Returns:

Name	Type	Description
p	ndarray of shape (n!, n)	Each row is a permutation.
s	ndarray of shape (n!,), optional	Signature (+1 or -1) of each permutation.

Source code in pyvoicebox/v_permutes.py

def v_permutes(n, return_sign=False) -> np.ndarray:
    """Generate all N! permutations of 0:N-1 in lexical order.

    Note: Returns 0-based indices (unlike MATLAB's 1-based).

    Parameters
    ----------
    n : int
        Number of elements to permute.
    return_sign : bool, optional
        If True, also return the signature of each permutation.

    Returns
    -------
    p : ndarray of shape (n!, n)
        Each row is a permutation.
    s : ndarray of shape (n!,), optional
        Signature (+1 or -1) of each permutation.
    """
    p = np.array([[0]])
    m = 1
    if n > 1:
        for a in range(1, n):
            nrows = (a + 1) * m
            q = np.zeros((nrows, a + 1), dtype=int)
            r = np.arange(1, a + 2)  # [1, 2, ..., a+1] then map to 0-based
            ix = np.arange(m)
            for b in range(a + 1):
                q[ix, 0] = b
                # Map p values through r (adjusting indices)
                q[ix, 1:a + 1] = r[p]
                r[b] = b
                ix = ix + m
            m = m * (a + 1)
            p = q

    if return_sign:
        s = 1 - 2 * (np.arange(1, m + 1) // 2 % 2)
        return p, s
    return p

v_sort

V_SORT - Sort with forward and inverse index.

v_sort

v_sort(
    a, descend=False, return_inverse=False
) -> tuple[ndarray, ndarray]

Sort array with forward and optional inverse index.

Parameters:

Name	Type	Description	Default
a	array_like	Input vector or matrix (sorted along axis 0 for matrices).	required
descend	bool	Sort in descending order.	False
return_inverse	bool	If True, also return the inverse index j such that a = b[j].	False

Returns:

Name	Type	Description
b	ndarray	Sorted array.
i	ndarray	Forward index: b = a[i].
j	ndarray (only if return_inverse=True)	Inverse index: a = b[j].

Source code in pyvoicebox/v_sort.py

def v_sort(a, descend=False, return_inverse=False) -> tuple[np.ndarray, np.ndarray]:
    """Sort array with forward and optional inverse index.

    Parameters
    ----------
    a : array_like
        Input vector or matrix (sorted along axis 0 for matrices).
    descend : bool, optional
        Sort in descending order.
    return_inverse : bool, optional
        If True, also return the inverse index j such that a = b[j].

    Returns
    -------
    b : ndarray
        Sorted array.
    i : ndarray
        Forward index: b = a[i].
    j : ndarray (only if return_inverse=True)
        Inverse index: a = b[j].
    """
    a = np.asarray(a)

    if a.ndim <= 1:
        if descend:
            i = np.argsort(-a)
        else:
            i = np.argsort(a)
        b = a[i]
        if return_inverse:
            j = np.empty_like(i)
            j[i] = np.arange(len(i))
            return b, i, j
        return b, i

    # Matrix case: sort each column
    if descend:
        i = np.argsort(-a, axis=0)
    else:
        i = np.argsort(a, axis=0)

    r, c = a.shape
    b = np.take_along_axis(a, i, axis=0)

    if return_inverse:
        j = np.zeros_like(i)
        for col in range(c):
            j[i[:, col], col] = np.arange(r)
        return b, i, j

    return b, i