Skip to content

utils module

This module contains some utility programs used by the package.

ScalarFunctionAndGradient: TypeAlias = Callable[[np.ndarray, list, Optional[bool]], float | tuple[float, np.ndarray]] module-attribute

Type of f(v, args, gr) that returns a scalar value and also a gradient if gr is True

bs_error_abort(msg='error, aborting')

Report error and exits with code 1

Parameters:

Name Type Description Default
msg str

specifies the error message

'error, aborting'

Returns:

Type Description
None

nothing

Source code in cupid_matching/utils.py
64
65
66
67
68
69
70
71
72
73
74
def bs_error_abort(msg: str = "error, aborting") -> None:
    """Report error and exits with code 1

    Args:
        msg: specifies the error message

    Returns:
        nothing
    """
    print_stars(f"{bs_name_func(3)}: {msg}")
    sys.exit(1)

bs_name_func(back=2)

Get the name of the current function, or further back in the stack

Parameters:

Name Type Description Default
back int

2 for the current function, 3 for the function that called it, etc

2

Returns:

Type Description
str

the name of the function requested

Source code in cupid_matching/utils.py
31
32
33
34
35
36
37
38
39
40
41
42
def bs_name_func(back: int = 2) -> str:
    """Get the name of the current function, or further back in the stack

    Args:
        back: 2 for the current function, 3 for the function that called it, etc

    Returns:
        the name of the function requested
    """
    stack = extract_stack()
    func_name: str = stack[-back][2]
    return func_name

check_gradient_scalar_function(fg, p, args, mode='central', EPS=1e-06)

Checks the gradient of a scalar function.

Parameters:

Name Type Description Default
fg Any

should return the scalar value, and the gradient if its gr argument is True

required
p np.ndarray

where we are checking the gradient

required
args list[Any]

other arguments passed to fg

required
mode str

"central" or "forward" derivatives

'central'
EPS float

the step for forward or central derivatives

1e-06

Returns:

Type Description
tuple[np.ndarray, np.ndarray]

the analytic and numeric gradients

Source code in cupid_matching/utils.py
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
def check_gradient_scalar_function(
    fg: Any,
    p: np.ndarray,
    args: list[Any],
    mode: str = "central",
    EPS: float = 1e-6,
) -> tuple[np.ndarray, np.ndarray]:
    """Checks the gradient of a scalar function.

    Args:
        fg: should return the scalar value, and the gradient if its `gr` argument is `True`
        p: where we are checking the gradient
        args: other arguments passed to `fg`
        mode: "central" or "forward" derivatives
        EPS: the step for forward or central derivatives

    Returns:
        the analytic and numeric gradients
    """
    # if not isinstance(fg, ScalarFunctionAndGradient):
    #     bs_error_abort("wrong type for fg")
    f0, f_grad = fg(p, args, True)

    print_stars("checking the gradient: analytic, numeric")

    g = np.zeros_like(p)
    sp: int = p.size
    if mode == "central":
        for i in range(sp):
            x = p[i]
            p1 = p.copy()
            p1[i] = x + EPS
            f_plus = fg(p1, args, False)
            p1[i] -= 2.0 * EPS
            f_minus = fg(p1, args, False)
            g[i] = (f_plus - f_minus) / (2.0 * EPS)
            print(f"{i}: {f_grad[i]}, {g[i]}")
    elif mode == "forward":
        for i in range(sp):
            x = p[i]
            p1 = p.copy()
            p1[i] = x + EPS
            f_plus = fg(p1, args, False)
            g[i] = (f_plus - f0) / EPS
            print(f"{i}: {f_grad[i]}, {g[i]}")
    else:
        bs_error_abort("mode must be 'central' or 'forward'")

    return f_grad, g

der_nppow(a, b)

evaluates the derivatives in a and b of element-by-element a^b

Parameters:

Name Type Description Default
a np.ndarray

input Numpy arrays

required
b int | float | np.ndarray

a Numpy array of the same shape, or a scalar

required

Returns:

Type Description
TwoArrays

a pair of two arrays of the same shape as a

Source code in cupid_matching/utils.py
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
def der_nppow(a: np.ndarray, b: int | float | np.ndarray) -> TwoArrays:
    """
    evaluates the derivatives in a and b of element-by-element $a^b$

    Args:
        a: input Numpy arrays
        b: a Numpy array of the same shape, or a scalar

    Returns:
        a pair of two arrays of the same shape as `a`
    """
    mina: float = np.min(a)
    if mina <= 0:
        print_stars("All elements of a must be positive in der_nppow!")
        sys.exit(1)

    if isinstance(b, (int, float)):
        a_pow_b = np.power(a, b)
        return (b * a_pow_b / a, a_pow_b * log(a))
    else:
        if a.shape != b.shape:
            print_stars("nppow: b is not a number or an array of the same shape as a!")
            sys.exit(1)
        avec = a.ravel()
        bvec = b.ravel()
        a_pow_b = avec**bvec
        der_wrt_a = a_pow_b * bvec / avec
        der_wrt_b = a_pow_b * nplog(avec)
        return der_wrt_a.reshape(a.shape), der_wrt_b.reshape(a.shape)

describe_array(v, name='The array')

Descriptive statistics on an array interpreted as a vector

Parameters:

Name Type Description Default
v np.ndarray

the array

required
name str

its name

'The array'

Returns:

Type Description
NamedTuple

a DescribeResult namedtuple

Source code in cupid_matching/utils.py
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def describe_array(v: np.ndarray, name: str = "The array") -> NamedTuple:
    """Descriptive statistics on an array interpreted as a vector

    Args:
        v: the array
        name: its name

    Returns:
        a `DescribeResult` namedtuple
    """
    print_stars(f"{name} has:")
    d = sts.describe(v, None)
    print(f"Number of elements: {d.nobs}")
    print(f"Minimum: {d.minmax[0]}")
    print(f"Maximum: {d.minmax[1]}")
    print(f"Mean: {d.mean}")
    print(f"Stderr: {sqrt(d.variance)}")
    return d

npexp(a, deriv=False, bigx=30.0, verbose=False)

C^2 extension of \exp(a) above bigx

Parameters:

Name Type Description Default
a np.ndarray

a Numpy array

required
deriv bool

if True, the first derivative is also returned

False
bigx float

an upper bound

30.0
verbose bool

whether diagnoses are printed

False

Returns:

Name Type Description
bigx np.ndarray | TwoArrays

upper bound \exp(a) C^2-extended above bigx,

np.ndarray | TwoArrays

with its derivative if deriv is True

Source code in cupid_matching/utils.py
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
def npexp(
    a: np.ndarray,
    deriv: bool = False,
    bigx: float = 30.0,
    verbose: bool = False,
) -> np.ndarray | TwoArrays:
    """
    $C^2$ extension of  $\\exp(a)$ above `bigx`

    Args:
        a: a Numpy array
        deriv: if `True`, the first derivative is also returned
        bigx: an upper bound
        verbose: whether diagnoses are printed

    Returns:
        bigx: upper bound $\\exp(a)$  $C^2$-extended above `bigx`,
        with its derivative if `deriv` is `True`
    """
    if np.max(a) < bigx:
        expa = np.exp(a)
        return (expa, expa) if deriv else expa
    else:
        exparr = np.exp(np.minimum(a, bigx))
        ebigx = exp(bigx)
        darr = a - bigx
        exparr_larger = ebigx * (1.0 + darr * (1.0 + 0.5 * darr))
        if verbose:
            n_large_args = np.sum(a > bigx)
            if n_large_args > 0:
                finals = "s" if n_large_args > 1 else ""
                print(
                    f"npexp: {n_large_args} argument{finals} larger than {bigx}: maxi = {np.max(a)}"
                )
        expa = np.where(a < bigx, exparr, exparr_larger)
        if deriv:
            der_exparr = np.exp(np.minimum(a, bigx))
            der_exparr_larger = ebigx * (1.0 + darr)
            der_expa = np.where(a < bigx, der_exparr, der_exparr_larger)
            return expa, der_expa
        else:
            return expa

nplog(a, deriv=False, eps=1e-30, verbose=False)

C^2 extension of \ln(a) below eps

Parameters:

Name Type Description Default
a np.ndarray

a Numpy array

required
deriv bool

if True, the first derivative is also returned

False
eps float

a lower bound

1e-30
verbose bool

whether diagnoses are printed

False

Returns:

Type Description
np.ndarray | tuple[np.ndarray, np.ndarray]

\ln(a) C^2-extended below eps,

np.ndarray | tuple[np.ndarray, np.ndarray]

with its derivative if deriv is True

Source code in cupid_matching/utils.py
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def nplog(
    a: np.ndarray,
    deriv: bool = False,
    eps: float = 1e-30,
    verbose: bool = False,
) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
    """$C^2$ extension of  $\\ln(a)$ below `eps`

    Args:
        a: a Numpy array
        deriv: if `True`, the first derivative is also returned
        eps: a lower bound
        verbose: whether diagnoses are printed

    Returns:
        $\\ln(a)$ $C^2$-extended below `eps`,
        with its derivative if `deriv` is `True`
    """
    if np.min(a) > eps:
        loga = np.log(a)
        return (loga, 1.0 / a) if deriv else loga
    else:
        logarreps = np.log(np.maximum(a, eps))
        logarr_smaller = log(eps) - (eps - a) * (3.0 * eps - a) / (2.0 * eps * eps)
        if verbose:
            n_small_args = np.sum(a < eps)
            if n_small_args > 0:
                finals = "s" if n_small_args > 1 else ""
                print(
                    f"nplog: {n_small_args} argument{finals} smaller than {eps}: mini = {np.min(a)}"
                )
        loga = np.where(a > eps, logarreps, logarr_smaller)
        if deriv:
            der_logarreps = 1.0 / np.maximum(a, eps)
            der_logarr_smaller = (2.0 * eps - a) / (eps * eps)
            der_loga = np.where(a > eps, der_logarreps, der_logarr_smaller)
            return loga, der_loga
        else:
            return loga

npmaxabs(a)

The maximum absolute value in an array

Parameters:

Name Type Description Default
a np.ndarray

the array

required

Returns:

Type Description
float

\max{\vert a \vert}

Source code in cupid_matching/utils.py
166
167
168
169
170
171
172
173
174
175
176
def npmaxabs(a: np.ndarray) -> float:
    """The maximum absolute value in an array

    Args:
        a: the array

    Returns:
        $\\max{\\vert a \\vert}$
    """
    maxi: float = np.max(np.abs(a))
    return maxi

nppow(a, b, deriv=False)

Evaluates a^b element-by-element

Parameters:

Name Type Description Default
a np.ndarray

a Numpy array

required
b int | float | np.ndarray

if an array, it should have the same shape as a

required
deriv bool

if True, the first derivatives wrt a and b are also returned

False

Returns:

Type Description
np.ndarray | ThreeArrays

an array of the same shape as a, and if deriv is True,

np.ndarray | ThreeArrays

the derivatives wrt a and b

Source code in cupid_matching/utils.py
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
def nppow(
    a: np.ndarray, b: int | float | np.ndarray, deriv: bool = False
) -> np.ndarray | ThreeArrays:
    """Evaluates $a^b$ element-by-element

    Args:
        a: a Numpy array
        b: if an array, it should have the same shape as `a`
        deriv: if `True`, the first derivatives wrt `a` and `b`
            are also returned

    Returns:
        an array of the same shape as `a`, and if `deriv` is `True`,
        the derivatives wrt `a` and `b`
    """
    mina = np.min(a)
    if mina <= 0:
        bs_error_abort("All elements of a must be positive!")

    if isinstance(b, (int, float)):
        a_pow_b = a**b
        if deriv:
            return (a**b, b * a_pow_b / a, a_pow_b * log(a))
        else:
            return a_pow_b
    else:
        if a.shape != b.shape:
            bs_error_abort(f"a has shape {a.shape} and b has shape {b.shape}")
        avec = a.ravel()
        bvec = b.ravel()
        a_pow_b = avec**bvec
        if deriv:
            der_wrt_a = a_pow_b * bvec / avec
            der_wrt_b = a_pow_b * nplog(avec)
            return (
                a_pow_b.reshape(a.shape),
                der_wrt_a.reshape(a.shape),
                der_wrt_b.reshape(a.shape),
            )
        else:
            return a_pow_b.reshape(a.shape)

nprepeat_col(v, n)

Creates a matrix with n columns, all equal to v

Parameters:

Name Type Description Default
v np.ndarray

a vector of size m

required
n int

the number of columns requested

required

:return: a matrix of shape (m, n)

Source code in cupid_matching/utils.py
137
138
139
140
141
142
143
144
145
146
147
148
def nprepeat_col(v: np.ndarray, n: int) -> np.ndarray:
    """Creates a matrix with `n` columns, all equal to `v`

    Args:
        v: a vector of size `m`
        n: the number of columns requested

    :return: a matrix of shape `(m, n)`
    """
    _ = test_vector(v, "nprepeat_col")
    w: np.ndarray = v[:, np.newaxis]
    return np.repeat(w, n, axis=1)

nprepeat_row(v, m)

Creates a matrix with m rows, all equal to v

Parameters:

Name Type Description Default
v np.ndarray

a vector of size n

required
m int

the number of rows requested

required

Returns:

Type Description
np.ndarray

a matrix of shape (m, n)

Source code in cupid_matching/utils.py
151
152
153
154
155
156
157
158
159
160
161
162
163
def nprepeat_row(v: np.ndarray, m: int) -> np.ndarray:
    """
    Creates a matrix with `m` rows, all equal to `v`

    Args:
        v: a vector of size `n`
        m: the number of rows requested

    Returns:
        a matrix of shape `(m, n)`
    """
    _ = test_vector(v, "nprepeat_row")
    return np.repeat(v[np.newaxis, :], m, axis=0)

print_stars(title=None, n=70)

Prints a starred line, or two around the title

Parameters:

Name Type Description Default
title Optional[str]

an optional title

None
n int

the number of stars on the line

70

Returns:

Type Description
None

nothing

Source code in cupid_matching/utils.py
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def print_stars(title: Optional[str] = None, n: int = 70) -> None:
    """Prints a starred line, or two around the title

    Args:
        title:  an optional title
        n: the number of stars on the line

    Returns:
        nothing
    """
    line_stars = "*" * n
    print()
    print(line_stars)
    if title:
        print(title.center(n))
        print(line_stars)
    print()

test_matrix(x, fun_name=None)

Tests that x is a matrix; aborts otherwise

Parameters:

Name Type Description Default
x Any

a potential matrix

required
fun_name Optional[str]

the name of the calling function

None

Returns:

Type Description
tuple[int, int]

the shape of x if it is a matrix

Source code in cupid_matching/utils.py
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def test_matrix(x: Any, fun_name: Optional[str] = None) -> tuple[int, int]:
    """Tests that `x` is a matrix; aborts otherwise

    Args:
        x: a potential matrix
        fun_name: the name of the calling function

    Returns:
        the shape of `x` if it is a matrix
    """
    fun_str: str = "" if fun_name is None else fun_name + ":"
    if not isinstance(x, np.ndarray):
        bs_error_abort(f"x in {fun_str} should be a Numpy array")
    ndims_x = x.ndim
    if ndims_x != 2:
        bs_error_abort(f"x in {fun_str} should have two dimensions, not {ndims_x}")
    shx: tuple[int, int] = x.shape
    return shx

test_vector(x, fun_name=None)

Tests that x is a vector; aborts otherwise

Parameters:

Name Type Description Default
x Any

a potential vector

required
fun_name Optional[str]

the name of the calling function

None

Returns:

Type Description
int

the size of x if it is a vector

Source code in cupid_matching/utils.py
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def test_vector(x: Any, fun_name: Optional[str] = None) -> int:
    """Tests that `x` is a vector; aborts otherwise

    Args:
        x: a potential vector
        fun_name: the name of the calling function

    Returns:
        the size of `x` if it is a vector
    """
    fun_str = ["" if fun_name is None else fun_name + ":"]
    if not isinstance(x, np.ndarray):
        bs_error_abort(f"{fun_str} x should be a Numpy array")
    ndims_x = x.ndim
    if ndims_x != 1:
        bs_error_abort(f"{fun_str} x should have one dimension, not {ndims_x}")
    sx: int = x.size
    return sx