Coverage for /opt/conda/lib/python3.14/site-packages/medil/sample.py: 95%
66 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-23 02:29 +0000
1"""Generate random minimum MeDIL causal model graph or parameters."""
3import numpy as np
4import numpy.typing as npt
5from numpy.random import default_rng
7from ._ecc_algorithms import _find_clique_min_cover
8from .models import GaussianMCM, NeuroCausalFactorAnalysis
11def mcm(
12 rng: np.random.Generator = default_rng(0),
13 parameterization: str = "Gaussian",
14 biadj: npt.NDArray = np.array([]),
15 **kwargs,
16) -> GaussianMCM | NeuroCausalFactorAnalysis:
17 """Randomly generate a minimum MeDIL causal model with parameters.
19 Parameters
20 ----------
21 rng : numpy.random.Generator, optional
22 Random number generator. Default is ``default_rng(0)``.
23 parameterization : str, optional
24 Either ``"Gaussian"`` (default) for a linear Gaussian model or
25 ``"VAE"`` for a randomly initialized masked VAE model.
26 biadj : ndarray, optional
27 Biadjacency matrix to use. If empty (default), one is generated
28 randomly using :func:`biadj` with any extra keyword arguments.
29 **kwargs
30 Additional keyword arguments passed to :func:`biadj` when
31 generating a random biadjacency matrix.
33 Returns
34 -------
35 GaussianMCM or NeuroCausalFactorAnalysis
36 A model with randomly generated structure and parameters, ready to
37 call :meth:`~medil.models.GaussianMCM.sample` on without fitting to data.
38 """
39 if biadj.size == 0:
40 biadj = _biadj(rng=rng, **kwargs)
41 if parameterization == "Gaussian":
42 mcm = GaussianMCM(biadj=biadj, rng=rng)
43 params = mcm.parameters
45 num_edges = biadj.sum()
46 weights = (rng.random(num_edges) * 1.5) + 0.5
47 weights[rng.choice((True, False), num_edges)] *= -1
48 params.biadj_weights = np.zeros_like(biadj, float)
49 params.biadj_weights[biadj] = weights
51 num_meas = biadj.shape[1]
52 params.error_means = rng.random(num_meas) * 2
53 params.error_means[rng.choice((True, False), num_meas)] *= -1
55 params.error_variances = (rng.random(num_meas) * 1.5) + 0.5
57 elif parameterization == "VAE":
58 try:
59 import torch
60 from ._vae import VariationalAutoencoder
61 except ImportError:
62 raise ImportError(
63 "parameterization='VAE' requires PyTorch. "
64 "Install it with: pip install medil[ncfa]"
65 )
66 model = NeuroCausalFactorAnalysis(biadj=biadj, rng=rng)
67 num_latent, num_meas = biadj.shape
68 biadj_tensor = torch.tensor(biadj.T, dtype=torch.float32)
69 vae = VariationalAutoencoder(
70 num_latent=num_latent,
71 num_meas=num_meas,
72 num_hidden_layers=model.hyperparams["num_hidden_layers"],
73 latent_width=model.hyperparams["latent_width"],
74 meas_width=model.hyperparams["meas_width"],
75 biadj=biadj_tensor,
76 encoder_hidden_dim=model.hyperparams["encoder_hidden_dim"],
77 ).to(model.device)
78 model.parameters.vae = vae
79 return model
81 else:
82 raise ValueError(f"Parameterization '{parameterization}' is invalid.")
84 return mcm
87def biadj(
88 num_meas: int,
89 density: float = 0.2,
90 one_pure_child: bool = True,
91 num_latent: int = 0,
92 rng: np.random.Generator = default_rng(0),
93) -> npt.NDArray:
94 """Randomly generate a biadjacency matrix for a minimum MeDIL causal model.
96 Parameters
97 ----------
98 num_meas : int
99 Number of measurement (observed) variables.
100 density : float, optional
101 Controls how many measurement variables share latent parents.
102 0 gives one latent per measurement (no sharing); 1 gives maximum
103 sharing. Default 0.2.
104 one_pure_child : bool, optional
105 If True (default), each latent variable has at least one measurement
106 variable that it is the sole parent of (the one-pure-child assumption).
107 If False, the graph is drawn from an Erdős–Rényi random graph over
108 observed variables and the minimum edge clique cover is computed.
109 num_latent : int, optional
110 Number of latent variables. Only used when ``one_pure_child=True``.
111 If 0 (default), drawn uniformly from ``[1, num_meas)``.
112 rng : numpy.random.Generator, optional
113 Random number generator. Default is ``default_rng(0)``.
115 Returns
116 -------
117 biadj : ndarray of shape (num_latent, num_meas), dtype bool
118 Boolean biadjacency matrix where ``biadj[i, j]`` is True iff
119 latent variable ``i`` is a parent of measurement variable ``j``.
120 """
121 if one_pure_child:
122 if num_latent == 0:
123 num_latent = rng.integers(1, num_meas)
124 if density is None:
125 density = rng.random()
127 # specify pure children/independent set
128 biadj = np.zeros((num_latent, num_meas), bool)
129 biadj[:, :num_latent] = np.eye(num_latent)
131 # every child gets a parent; specifically L_0, until the
132 # within-column perm below using np.permuted
133 biadj[0, num_latent:] = True
135 # randomly fill in remaining density * (num_meas - num_latent)
136 # * (num_latent - 1) edges
137 max_num_edges = (num_meas - num_latent) * (num_latent - 1)
138 num_edges = np.round(max_num_edges * density).astype(int)
140 edges = np.zeros(max_num_edges, bool)
141 edges[:num_edges] = True
142 edges = rng.permutation(edges).reshape(num_latent - 1, num_meas - num_latent)
144 biadj[1:][:, num_latent:] = edges
146 nonpure_children = biadj[:, num_latent:]
147 biadj[:, num_latent:] = rng.permuted(nonpure_children, axis=0)
149 # change child order, so pure children aren't first
150 biadj = rng.permutation(biadj, axis=1)
152 else:
153 if num_latent != 0:
154 msg = "`num_latent` can only be specified when `one_pure_child==True`."
155 raise ValueError(msg)
157 udg = np.zeros((num_meas, num_meas), bool)
159 max_edges = (num_meas * (num_meas - 1)) // 2
160 num_edges = np.round(density * max_edges).astype(int)
162 edges = np.ones(max_edges)
163 edges[num_edges:] = 0
165 udg[np.triu_indices(num_meas, k=1)] = rng.permutation(edges)
166 udg += udg.T
167 np.fill_diagonal(udg, True)
169 # find latent connections (minimum edge clique cover)
170 biadj = _find_clique_min_cover(udg).astype(bool)
172 return biadj
175_biadj = biadj