API Reference

Auto-generated from the source docstrings. For a task-oriented walkthrough with examples, see the API Guide.

All symbols below are importable directly from the top-level logit_graph package, e.g. from logit_graph import simulate_graph.

The model — Logistic Random Graph (LG)

There is a single model, the Logistic Random Graph (LG) of Ottoni, Takahashi & Fujita (2026); the entries below are its computational forms and the general-pairwise-features extension (§3.7). All score an edge with a logistic link and differ only in what enters the linear predictor η_ij.

Form	Edge log-odds η_ij	Packaged as
Equilibrium sampler	σ + α·[g(S_i)+g(S_j)] at stationarity (Gibbs)	simulate_graph, GraphModel
Iterative generation	same, with S evaluated on the lagged graph G(t-1) (Eq. 3.5)	LG generation & estimation
General pairwise features (§3.7)	σ + Σₖ θₖ·φₖ(i,j) — degree term plus extra features	Unified LG

• σ — baseline connection propensity (baseline log-odds; typically negative).
• α ≥ 0 — strength of the neighborhood degree effect on g(S_i)+g(S_j), g(s)=log(1+s), S_i = Σ_{l∈N_d(i)} k_l the neighborhood degree sum over the d-hop ball.
• d — the neighborhood radius (hops). d=0 uses each node's own degree only (S_i = k_i); d=1 adds the neighbors' degrees; larger d reaches further. Selected from data by AIC.
• φₖ (extra features) — appended columns: exogenous / pair-local ones (e.g. block membership via community_feature) are directly estimable; globally-coupled ones (e.g. latent proximity via latent_feature, or the degree term itself) are made consistent by the predetermined-predictor (temporal) form that reads them off G(t-1). The unified five-feature LG model of the thesis is the instance with degree + coarse/fine community (γc, γf) + latent proximity (λ).

Paper-consistent sampler & estimator

simulate_graph

simulate_graph(
    n: int,
    d: int,
    sigma: Optional[float] = None,
    *,
    beta: Optional[float] = None,
    n_iter: int,
    feature_mode: FeatureMode = "incremental",
    target_density: float = 0.1,
    signal: float = 0.5,
    seed: Optional[int] = None,
    return_meta: bool = False,
    collect_feature_mean: bool = False,
    materialize_adjacency: bool = True,
    feature_mode_est: Optional[FeatureMode] = None,
    adaptive_stopping: bool = False,
    adaptive_check_interval: int = 20000,
    adaptive_patience: int = 3,
    adaptive_cv_tol: float = 0.02,
    adaptive_min_iter: int = 20000
) -> np.ndarray | tuple[np.ndarray, dict[str, float]]

Generate a graph at fixed sigma. d=0 uses a direct ER sample at p=expit(sigma) (exact equilibrium, no degree feedback); d>=1 runs leave-one-out Gibbs with the feature coefficient fixed at beta=1 (matching the paper offset estimator).

Source code in src/logit_graph/experiments/sweeps.py

def simulate_graph(
    n: int,
    d: int,
    sigma: Optional[float] = None,
    *,
    beta: Optional[float] = None,
    n_iter: int,
    feature_mode: FeatureMode = "incremental",
    target_density: float = 0.10,
    signal: float = 0.5,
    seed: Optional[int] = None,
    return_meta: bool = False,
    collect_feature_mean: bool = False,
    materialize_adjacency: bool = True,
    feature_mode_est: Optional[FeatureMode] = None,
    adaptive_stopping: bool = False,
    adaptive_check_interval: int = 20_000,
    adaptive_patience: int = 3,
    adaptive_cv_tol: float = 0.02,
    adaptive_min_iter: int = 20_000,
) -> np.ndarray | tuple[np.ndarray, dict[str, float]]:
    """Generate a graph at fixed sigma. d=0 uses a direct ER sample at p=expit(sigma)
    (exact equilibrium, no degree feedback); d>=1 runs leave-one-out Gibbs with the feature
    coefficient fixed at beta=1 (matching the paper offset estimator)."""
    if d == 0:
        if sigma is None:
            sigma = math.log(target_density / (1.0 - target_density))
        adj = _sample_er_at_sigma(n, sigma, seed)
        if not return_meta:
            return adj
        meta = {
            "sigma": float(sigma),
            "beta": 0.0,
            "density": _graph_density(adj),
            "feature_mean": 0.0,
        }
        return adj, meta

    if seed is not None and beta is None:
        if sigma is not None:
            beta = _calibrate_beta_given_sigma(
                n, d, sigma, target_density, signal, feature_mode, seed,
            )
        else:
            sigma_cal, beta_cal = _calibrate_beta(
                n, d, target_density, signal, feature_mode, seed,
            )
            sigma = sigma_cal
            beta = beta_cal

    assert sigma is not None
    use_beta = beta if beta is not None else 1.0
    gen_mode: FeatureMode = feature_mode

    from scipy.special import expit as _expit
    # Warm-start near equilibrium: use expit(sigma) for very sparse cases,
    # but never below a minimum that lets d>=1 features (ball overlaps) be
    # populated enough to inform Gibbs updates.
    er_init = float(np.clip(_expit(sigma), 0.02, 0.5))
    n_iter_used = n_iter

    gm = graph_mod.GraphModel(
        n=n,
        d=d,
        sigma=sigma,
        alpha=1.0,
        beta=float(use_beta),
        er_p=er_init,
        layer2=True,
        feature_mode=gen_mode,
        seed=seed,
    )

    if adaptive_stopping and d >= 1:
        from ..lg_features_fast import FastGibbsGraph

        fg = FastGibbsGraph(
            n,
            d,
            sigma,
            er_p=er_init,
            rng=gm._rng,
            feature_mode=gen_mode,
            alpha=1.0,
            beta=float(use_beta),
            adj=gm.graph,
        )
        csr_rows, n_iter_used = _run_adaptive_gibbs(
            fg,
            gm._rng,
            n_iter,
            check_interval=adaptive_check_interval,
            patience=adaptive_patience,
            cv_tol=adaptive_cv_tol,
            min_iter=adaptive_min_iter,
        )
    else:
        gm.populate_edges_baseline(
            warm_up=0,
            max_iterations=n_iter,
            patience=10,
            check_interval=10**9,
            fast_mode=True,
        )
        csr_rows = getattr(gm, "_csr_rows", None)
    from ..lg_features_fast import adj_from_rows, density_from_rows

    if csr_rows is not None:
        adj = adj_from_rows(csr_rows, n) if materialize_adjacency else None
    else:
        adj = gm.graph.copy()

    if not return_meta:
        if adj is None:
            adj = adj_from_rows(csr_rows, n)
        return adj

    density = (
        density_from_rows(csr_rows, n) if csr_rows is not None else _graph_density(adj)
    )
    if collect_feature_mean:
        feat_seed = seed if seed is not None else 0
        if csr_rows is not None and d >= 1:
            from ..lg_features_fast import build_pair_dataset_from_rows

            offsets, _ = build_pair_dataset_from_rows(
                csr_rows, d, mode=gen_mode,
            )
            feat_mean = float(np.mean(offsets))
        else:
            feat_mean = _mean_pair_feature(adj, d, gen_mode, feat_seed)
    else:
        feat_mean = 0.0

    meta = {
        "sigma": float(sigma),
        "beta": float(use_beta),
        "density": density,
        "feature_mean": feat_mean,
        "n_iter_used": float(n_iter_used),
    }
    if csr_rows is not None:
        meta["csr_rows"] = csr_rows
    return adj, meta

select_d_ensemble

select_d_ensemble(
    graphs: list[ndarray],
    d_candidates: list[int],
    feature_mode: FeatureMode,
    extra_penalty_per_d: float = 0.0,
    *,
    csr_rows_list: Optional[list[Optional[list]]] = None
) -> tuple[int, dict[int, dict[str, float]]]

Source code in src/logit_graph/experiments/sweeps.py

def select_d_ensemble(
    graphs: list[np.ndarray],
    d_candidates: list[int],
    feature_mode: FeatureMode,
    extra_penalty_per_d: float = 0.0,
    *,
    csr_rows_list: Optional[list[Optional[list]]] = None,
) -> tuple[int, dict[int, dict[str, float]]]:
    from ..lg_features_fast import build_multi_d_pair_datasets_fast

    d_sorted = sorted(d_candidates)
    offsets_acc: dict[int, list[np.ndarray]] = {d: [] for d in d_candidates}
    labels_acc: list[np.ndarray] = []

    for idx, g in enumerate(graphs):
        rows = None if csr_rows_list is None else csr_rows_list[idx]
        labels, offsets_by_d = build_multi_d_pair_datasets_fast(
            g, d_sorted, mode=feature_mode, rows=rows,
        )
        labels_acc.append(labels)
        for d in d_candidates:
            offsets_acc[d].append(offsets_by_d[d])

    labels_all = np.concatenate(labels_acc)
    stats = {
        d: _ensemble_aic_stats(
            np.concatenate(offsets_acc[d]),
            labels_all,
            d,
            extra_penalty=extra_penalty_per_d * d,
        )
        for d in d_candidates
    }
    valid = {d: s for d, s in stats.items() if np.isfinite(s["aic"])}
    best = min(valid, key=lambda d: valid[d]["aic"]) if valid else d_candidates[0]
    return best, stats

estimate_sigma_from_graph

estimate_sigma_from_graph(
    adj: Optional[ndarray],
    d: int,
    feature_mode: FeatureMode = "incremental",
    beta: float = 1.0,
    *,
    csr_rows: Optional[list] = None
) -> float

Source code in src/logit_graph/experiments/sweeps.py

def estimate_sigma_from_graph(
    adj: Optional[np.ndarray],
    d: int,
    feature_mode: FeatureMode = "incremental",
    beta: float = 1.0,
    *,
    csr_rows: Optional[list] = None,
) -> float:
    del beta
    if csr_rows is not None and d >= 1:
        from ..lg_features_fast import build_pair_dataset_from_rows
        from ..offset_logit import fit_offset_logit_fast

        offsets, labels = build_pair_dataset_from_rows(csr_rows, d, mode=feature_mode)
        sigma_hat, _ = fit_offset_logit_fast(offsets, labels)
        return float(sigma_hat)
    if adj is None:
        raise ValueError("adj required when csr_rows is not provided")
    est = LogitRegEstimator(adj, d=d, layer2=True, feature_mode=feature_mode)
    stats = est.compute_aic(d_est=d, feature_mode=feature_mode)
    return float(stats["sigma_hat"])

High-level API (scikit-learn style)

LogitGraphFitter

LogitGraphFitter(
    d: int = 0,
    n_iteration: int = 10000,
    warm_up: int = 500,
    patience: int = 2000,
    dist_type: str = "KL",
    edge_delta: Optional[float] = None,
    min_gic_threshold: float = 5,
    check_interval: int = 50,
    verbose: bool = True,
    er_p: float = 0.05,
    init_graph: Optional[Graph] = None,
)

Fits a single Logit Graph model to a real graph, following a scikit-learn-like API.

Initialize the LogitGraphFitter: d (latent depth), n_iteration cap, warm_up, patience (checks without improvement before stopping), dist_type for GIC, edge_delta / min_gic_threshold convergence thresholds, check_interval, verbose.

Source code in src/logit_graph/simulation.py

def __init__(
    self,
    d: int = 0,
    n_iteration: int = 10000,
    warm_up: int = 500,
    patience: int = 2000,
    dist_type: str = 'KL',
    edge_delta: Optional[float] = None,
    min_gic_threshold: float = 5,
    check_interval: int = 50,
    verbose: bool = True,
    er_p: float = 0.05,
    init_graph: Optional[nx.Graph] = None,
) -> None:
    """Initialize the LogitGraphFitter: d (latent depth), n_iteration cap, warm_up,
    patience (checks without improvement before stopping), dist_type for GIC,
    edge_delta / min_gic_threshold convergence thresholds, check_interval, verbose."""
    self.d = d
    self.n_iteration = n_iteration
    self.warm_up = warm_up
    self.er_p = er_p
    self.patience = patience
    self.dist_type = dist_type
    self.edge_delta = edge_delta
    self.min_gic_threshold = min_gic_threshold
    self.check_interval = check_interval
    self.verbose = verbose
    self.init_graph = init_graph

    self.fitted_graph = None
    self.metadata = {}

fit

fit(original_graph: Graph) -> LogitGraphFitter

Fit the Logit Graph model to original_graph; returns self with the fitted_graph and metadata attributes populated.

Source code in src/logit_graph/simulation.py

def fit(self, original_graph: nx.Graph) -> LogitGraphFitter:
    """Fit the Logit Graph model to ``original_graph``; returns self with the
    ``fitted_graph`` and ``metadata`` attributes populated."""
    # Ensure we work with an undirected view for consistent edge counting and spectra
    # Many datasets load as directed; our fitter/populator assumes undirected adjacency.
    undirected_graph = original_graph.to_undirected()

    if self.verbose:
        print(f"\n{'='*20} Processing Graph {'='*20}")
        print(f"Original graph - Nodes: {undirected_graph.number_of_nodes()}, Edges: {undirected_graph.number_of_edges()}")

    self.metadata = {
        'original_nodes': undirected_graph.number_of_nodes(),
        'original_edges': undirected_graph.number_of_edges(),
        'fit_success': False,
        'error_message': None,
    }

    try:
        # Build adjacency from the undirected graph to avoid double-counting edges
        adj_matrix = nx.to_numpy_array(undirected_graph)

        best_graph_arr, sigma, gic_val, spectrum_diffs, edge_diffs, best_iter, all_graphs, gic_values = self._generate_graph(adj_matrix)

        self.fitted_graph = nx.from_numpy_array(best_graph_arr)

        self.metadata.update({
            'fit_success': True,
            'sigma': sigma,
            'gic_value': gic_val,
            'best_iteration': best_iter,
            'fitted_nodes': self.fitted_graph.number_of_nodes(),
            'fitted_edges': self.fitted_graph.number_of_edges(),
            'spectrum_diffs': spectrum_diffs,
            'edge_diffs': edge_diffs,
            'gic_values': gic_values,
        })

        if self.verbose:
            print(f"Fitting successful - GIC: {self.metadata['gic_value']:.4f}, Best iteration: {self.metadata['best_iteration']}")
            print(f"Fitted graph - Nodes: {self.metadata['fitted_nodes']}, Edges: {self.metadata['fitted_edges']}")

        del all_graphs
        gc.collect()

    except Exception as e:
        print(f"Error fitting graph: {e}")
        self.metadata['error_message'] = str(e)
        self.fitted_graph = None

    return self

LogitGraphSimulation

LogitGraphSimulation(
    n: int,
    d: int,
    sigma: float,
    alpha: float = 1.0,
    beta: float = 1.0,
    er_p: float = 0.05,
    n_iteration: int = 10000,
    warm_up: int = 500,
    patience: int = 2000,
    check_interval: int = 50,
    edge_cv_tol: float = 0.02,
    spectrum_cv_tol: float = 0.02,
    verbose: bool = True,
    init_graph: Optional[Graph] = None,
    layer2: bool = True,
    feature_mode: FeatureMode = "bounded",
    fast_mode: bool = False,
)

Simulate a Logit Graph directly from provided parameters (scikit-learn-like API): sim = LogitGraphSimulation(n=100, d=2, sigma=1.0, alpha=1.0, beta=1.0); sim.simulate(), then read sim.simulated_graph and sim.metadata.

Initialize the LogitGraphSimulation: model params (n, d, sigma, alpha, beta, er_p) and run params (n_iteration, warm_up, patience window, check_interval, edge/spectrum CV tolerances, verbose, init_graph).

Source code in src/logit_graph/simulation.py

def __init__(
    self,
    n: int,
    d: int,
    sigma: float,
    alpha: float = 1.0,
    beta: float = 1.0,
    er_p: float = 0.05,
    n_iteration: int = 10000,
    warm_up: int = 500,
    patience: int = 2000,
    check_interval: int = 50,
    edge_cv_tol: float = 0.02,
    spectrum_cv_tol: float = 0.02,
    verbose: bool = True,
    init_graph: Optional[nx.Graph] = None,
    layer2: bool = True,
    feature_mode: FeatureMode = "bounded",
    fast_mode: bool = False,
) -> None:
    """Initialize the LogitGraphSimulation: model params (n, d, sigma, alpha, beta,
    er_p) and run params (n_iteration, warm_up, patience window, check_interval,
    edge/spectrum CV tolerances, verbose, init_graph)."""
    self.n = n
    self.d = d
    self.sigma = sigma
    self.alpha = alpha
    self.beta = beta
    self.er_p = er_p
    self.n_iteration = n_iteration
    self.warm_up = warm_up
    self.patience = patience
    self.check_interval = check_interval
    self.edge_cv_tol = edge_cv_tol
    self.spectrum_cv_tol = spectrum_cv_tol
    self.verbose = verbose
    self.init_graph = init_graph
    self.layer2 = layer2
    self.feature_mode = feature_mode
    self.fast_mode = fast_mode

    self.simulated_graph = None
    self.metadata = {}

simulate

simulate() -> LogitGraphSimulation

Run the simulation, storing the graph and metadata. d=0 samples directly from ER(p=expit(sigma)) (exact equilibrium); d>=1 warm-starts the Gibbs chain near expit(sigma) clipped to [0.02, 0.5], unless init_graph or a non-default er_p overrides.

Source code in src/logit_graph/simulation.py

def simulate(self) -> LogitGraphSimulation:
    """Run the simulation, storing the graph and metadata. d=0 samples directly from
    ER(p=expit(sigma)) (exact equilibrium); d>=1 warm-starts the Gibbs chain near
    expit(sigma) clipped to [0.02, 0.5], unless init_graph or a non-default er_p overrides."""
    if self.verbose:
        print(f"\n{'='*20} Simulating Logit Graph {'='*20}")
        print(f"Parameters - n: {self.n}, d: {self.d}, sigma: {self.sigma:.4f}, alpha: {self.alpha:.4f}, beta: {self.beta:.4f}, er_p: {self.er_p}")

    self.metadata = {
        'simulate_success': False,
        'error_message': None,
        'n': self.n,
        'd': self.d,
        'sigma': float(self.sigma),
        'alpha': float(self.alpha),
        'beta': float(self.beta),
        'er_p': float(self.er_p),
    }

    try:
        if self.d == 0:
            final_graph_arr = _direct_er_at_sigma(self.n, self.sigma, seed=None)
            self.simulated_graph = nx.from_numpy_array(final_graph_arr)
            self.metadata.update({
                'simulate_success': True,
                'iterations_ran': 0,
                'final_nodes': self.simulated_graph.number_of_nodes(),
                'final_edges': self.simulated_graph.number_of_edges(),
                'final_spectrum': None,
                'sampler': 'direct_er',
            })
            if self.verbose:
                print(
                    f"d=0: direct ER at p=expit({self.sigma:.3f})={expit(self.sigma):.4f} — "
                    f"Nodes={self.metadata['final_nodes']}, Edges={self.metadata['final_edges']}"
                )
            return self

        user_override = self.init_graph is not None or not math.isclose(self.er_p, 0.05)
        warm_start_p = self.er_p if user_override else _warm_start_er_p(self.sigma)

        graph_model = graph.GraphModel(
            n=self.n, d=self.d, sigma=self.sigma, alpha=self.alpha, beta=self.beta,
            er_p=warm_start_p, init_graph=self.init_graph,
            layer2=self.layer2, feature_mode=self.feature_mode,
        )

        graphs, spectra = graph_model.populate_edges_baseline(
            warm_up=self.warm_up, max_iterations=self.n_iteration,
            patience=self.patience, check_interval=self.check_interval,
            edge_cv_tol=self.edge_cv_tol, spectrum_cv_tol=self.spectrum_cv_tol,
            fast_mode=self.fast_mode,
        )

        final_graph_arr = graphs[-1] if graphs else graph_model.graph
        self.simulated_graph = nx.from_numpy_array(final_graph_arr)

        self.metadata.update({
            'simulate_success': True,
            'iterations_ran': max(0, len(graphs) - 1),
            'final_nodes': self.simulated_graph.number_of_nodes(),
            'final_edges': self.simulated_graph.number_of_edges(),
            'final_spectrum': spectra,
            'warm_start_er_p': float(warm_start_p),
            'sampler': 'gibbs_layer2',
        })

        if self.verbose:
            print(
                f"Simulation successful (warm_start_p={warm_start_p:.4f}) — "
                f"Nodes: {self.metadata['final_nodes']}, Edges: {self.metadata['final_edges']}"
            )

    except Exception as e:
        print(f"Error simulating logit graph: {e}")
        self.metadata['error_message'] = str(e)
        self.simulated_graph = None

    return self

GraphModelComparator

GraphModelComparator(
    d_list: list[int],
    lg_params: dict[str, Any],
    other_model_n_runs: int = 2,
    other_model_params: Optional[
        list[dict[str, Any]]
    ] = None,
    dist_type: str = "KL",
    verbose: bool = True,
    other_models: Optional[list[str]] = None,
    other_model_grid_points: int = 5,
    random_state: Optional[int] = 42,
)

Compares Logit Graph with other random graph models, with a scikit-learn-like API.

Initialize the GraphModelComparator: d_list to test, lg_params (LG fitting), other_model_n_runs / other_model_params for the baselines, dist_type for GIC, verbose, and random_state (None for non-reproducible runs).

Source code in src/logit_graph/simulation.py

def __init__(
    self,
    d_list: list[int],
    lg_params: dict[str, Any],
    other_model_n_runs: int = 2,
    other_model_params: Optional[list[dict[str, Any]]] = None,
    dist_type: str = 'KL',
    verbose: bool = True,
    other_models: Optional[list[str]] = None,
    other_model_grid_points: int = 5,
    random_state: Optional[int] = 42,
) -> None:
    """Initialize the GraphModelComparator: d_list to test, lg_params (LG fitting),
    other_model_n_runs / other_model_params for the baselines, dist_type for GIC,
    verbose, and random_state (None for non-reproducible runs)."""
    self.d_list = d_list
    self.lg_params = lg_params
    self.other_model_n_runs = other_model_n_runs
    self.random_state = random_state
    # ``other_model_params`` may be None (per-model defaults looked up by name in
    # ``_fit_other_models``; recommended), a list (positional pairing with
    # ``other_models``; legacy), or a dict keyed by model name (explicit overrides).
    self.other_model_params = other_model_params
    self.dist_type = dist_type
    self.verbose = verbose
    # Allow selecting which other models to evaluate and how dense the parameter grid should be
    # Default excludes GRG for speed; user can include it by passing other_models.
    self.other_models = other_models if other_models is not None else ["ER", "WS", "BA"]
    self.other_model_grid_points = other_model_grid_points

    self.summary_df = None
    self.fitted_graphs_data = {}

compare

compare(
    original_graph: Graph, graph_filepath: str
) -> GraphModelComparator

Fit LG and the other models to original_graph (graph_filepath used only for logging) and compare them; returns self with summary_df and fitted_graphs_data populated.

Source code in src/logit_graph/simulation.py

def compare(self, original_graph: nx.Graph, graph_filepath: str) -> GraphModelComparator:
    """Fit LG and the other models to ``original_graph`` (``graph_filepath`` used only
    for logging) and compare them; returns self with ``summary_df`` and
    ``fitted_graphs_data`` populated."""
    if self.verbose:
        print(f"\n{'='*30} Processing Graph: {os.path.basename(graph_filepath)} {'='*30}")

    if self.random_state is not None:
        random.seed(self.random_state)
        np.random.seed(self.random_state)

    self.fitted_graphs_data = {
        'Original': {
            'graph': original_graph,
            'metadata': {'fit_success': True, 'param': 'N/A', 'gic_value': np.nan}
        }
    }
    adj_matrix = nx.to_numpy_array(original_graph)

    # 1. Fit Logit Graph (LG) model, finding the best `d`
    self._fit_best_lg(adj_matrix)

    # 2. Fit other random graph models
    self._fit_other_models(original_graph)

    # 3. Calculate attributes and build the summary DataFrame
    self._build_summary_df(graph_filepath)

    return self

estimate_sigma_only

estimate_sigma_only(
    graph_input: Union[ndarray, Graph],
    d: int,
    max_pairs: Optional[int] = None,
    feature_mode: FeatureMode = "incremental",
    random_state: Optional[int] = None,
    verbose: bool = False,
    max_edges: Optional[int] = None,
    max_non_edges: Optional[int] = None,
    l1_wt: float = 1,
    alpha: float = 0,
) -> tuple[float, Any]

Estimate sigma via the leave-one-out (full conditional) offset logit (paper formulation), using :func:build_pair_dataset (layer2=True) so it stays consistent with generation. Returns (sigma_hat, fit_result); legacy max_edges/max_non_edges/l1_wt/alpha unused.

Source code in src/logit_graph/simulation.py

def estimate_sigma_only(
    graph_input: Union[np.ndarray, nx.Graph],
    d: int,
    max_pairs: Optional[int] = None,
    feature_mode: FeatureMode = "incremental",
    random_state: Optional[int] = None,
    verbose: bool = False,
    # Deprecated kwargs kept for backwards compatibility ------------------
    max_edges: Optional[int] = None,
    max_non_edges: Optional[int] = None,
    l1_wt: float = 1,
    alpha: float = 0,
) -> tuple[float, Any]:
    """Estimate sigma via the leave-one-out (full conditional) offset logit (paper formulation),
    using :func:`build_pair_dataset` (layer2=True) so it stays consistent with generation.
    Returns ``(sigma_hat, fit_result)``; legacy max_edges/max_non_edges/l1_wt/alpha unused."""
    del l1_wt, alpha  # accepted for backwards compatibility, unused

    if max_pairs is None and (max_edges is not None or max_non_edges is not None):
        max_pairs = (max_edges or 0) + (max_non_edges or 0) or None

    adj = _as_adj(graph_input)

    est = estimator.LogitRegEstimator(
        adj, d=d, layer2=True, feature_mode=feature_mode, verbose=verbose,
    )
    offsets, labels_arr = build_pair_dataset(
        adj, d=d, mode=feature_mode, layer2=True,
        max_pairs=max_pairs, seed=random_state,
    )
    result = est._fit_offset_logit(offsets, np.asarray(labels_arr, dtype=int))
    sigma = float(result.params[0])
    return sigma, result

estimate_sigma_many

estimate_sigma_many(
    graph_input: Union[ndarray, Graph],
    d: int,
    n_repeats: int = 30,
    max_pairs: Optional[int] = None,
    feature_mode: FeatureMode = "incremental",
    seed: int = 42,
    verbose: bool = False,
    max_edges: Optional[int] = None,
    max_non_edges: Optional[int] = None,
    l1_wt: float = 1,
    alpha: float = 0,
) -> list[float]

Repeat the leave-one-out sigma estimation n_repeats times with different seeds, returning a list of sigma_hat. Variability comes from random pair sampling when max_pairs is set; with max_pairs=None on a small graph every repetition is identical.

Source code in src/logit_graph/simulation.py

def estimate_sigma_many(
    graph_input: Union[np.ndarray, nx.Graph],
    d: int,
    n_repeats: int = 30,
    max_pairs: Optional[int] = None,
    feature_mode: FeatureMode = "incremental",
    seed: int = 42,
    verbose: bool = False,
    # Deprecated kwargs kept for backwards compatibility ------------------
    max_edges: Optional[int] = None,
    max_non_edges: Optional[int] = None,
    l1_wt: float = 1,
    alpha: float = 0,
) -> list[float]:
    """Repeat the leave-one-out sigma estimation ``n_repeats`` times with different seeds, returning
    a list of ``sigma_hat``. Variability comes from random pair sampling when ``max_pairs``
    is set; with ``max_pairs=None`` on a small graph every repetition is identical."""
    del l1_wt, alpha
    if max_pairs is None and (max_edges is not None or max_non_edges is not None):
        max_pairs = (max_edges or 0) + (max_non_edges or 0) or None

    sigmas: list[float] = []
    for r in tqdm(range(int(n_repeats))):
        rs = None if seed is None else (seed + r)
        sigma_r, _result = estimate_sigma_only(
            graph_input, d=d, max_pairs=max_pairs,
            feature_mode=feature_mode, random_state=rs, verbose=verbose,
        )
        sigmas.append(sigma_r)
    return sigmas

calculate_graph_attributes

calculate_graph_attributes(
    graph_to_analyze: Optional[Graph],
) -> dict[str, Any]

Calculate various graph attributes for a given graph.

Source code in src/logit_graph/simulation.py

def calculate_graph_attributes(graph_to_analyze: Optional[nx.Graph]) -> dict[str, Any]:
    """Calculate various graph attributes for a given graph."""
    if graph_to_analyze is None or graph_to_analyze.number_of_nodes() == 0:
        return {attr: np.nan for attr in [
            'nodes', 'edges', 'density', 'avg_clustering', 'avg_path_length', 'diameter',
            'assortativity', 'num_components', 'largest_component_size'
        ]}

    attrs = {'nodes': graph_to_analyze.number_of_nodes(), 'edges': graph_to_analyze.number_of_edges()}

    try:
        attrs['density'] = nx.density(graph_to_analyze)
        attrs['avg_clustering'] = nx.average_clustering(graph_to_analyze)
        attrs['assortativity'] = nx.degree_assortativity_coefficient(graph_to_analyze)

        components = list(nx.connected_components(graph_to_analyze))
        attrs['num_components'] = len(components)
        if components:
            largest_cc = max(components, key=len)
            attrs['largest_component_size'] = len(largest_cc)
            subgraph = graph_to_analyze.subgraph(largest_cc)
            if len(largest_cc) > 1:
                attrs['avg_path_length'] = nx.average_shortest_path_length(subgraph)
                attrs['diameter'] = nx.diameter(subgraph)
            else:
                attrs['avg_path_length'] = 0
                attrs['diameter'] = 0
        else:
             attrs['largest_component_size'] = 0
             attrs['avg_path_length'] = np.nan
             attrs['diameter'] = np.nan

    except Exception as e:
        print(f"Could not calculate all attributes: {e}")

    return attrs

Core model & estimator

GraphModel

GraphModel(
    n: int,
    d: int,
    sigma: float,
    alpha: float = 1,
    beta: float = 1,
    er_p: float = 0.05,
    init_graph: Optional[Graph] = None,
    layer2: bool = True,
    feature_mode: FeatureMode = "incremental",
    seed: Optional[int] = None,
)

Source code in src/logit_graph/graph.py

def __init__(
    self,
    n: int,
    d: int,
    sigma: float,
    alpha: float = 1,
    beta: float = 1,
    er_p: float = 0.05,
    init_graph: Optional[nx.Graph] = None,
    layer2: bool = True,
    feature_mode: FeatureMode = "incremental",
    seed: Optional[int] = None,
) -> None:
    # NOTE: feature_mode defaults to "incremental" so d=0 has zero pair feature
    # (pure Erdős–Rényi at p=expit(sigma)) and d>=1 uses a bounded, identifying
    # feature; the legacy "bounded" mode adds degree feedback even at d=0 and saturates.
    self.n = n
    self.d = d
    self.sigma = sigma
    self.alpha = alpha
    self.beta = beta
    self.er_p = er_p
    self.layer2 = layer2
    self.feature_mode = feature_mode
    self._adj_only: bool = False

    self._rng = np.random.default_rng(seed)

    if init_graph is not None and isinstance(init_graph, nx.Graph):
        self.graph = nx.to_numpy_array(init_graph)
    else:
        self.graph = self.generate_small_er_graph(n, p=er_p)

    # Cached state — kept in sync by add_remove_edge
    self._init_cache()

generate_small_er_graph

generate_small_er_graph(n: int, p: float) -> np.ndarray

Vectorized ER sample (symmetric, uses self._rng).

Source code in src/logit_graph/graph.py

def generate_small_er_graph(self, n: int, p: float) -> np.ndarray:
    """Vectorized ER sample (symmetric, uses ``self._rng``)."""
    upper = self._rng.random((n, n)) < p
    upper = np.triu(upper, k=1)
    return (upper | upper.T).astype(float)

calculate_spectrum classmethod

calculate_spectrum(graph: ndarray) -> np.ndarray

Sorted eigenvalues of the Laplacian L = D - A.

Source code in src/logit_graph/graph.py

@classmethod
def calculate_spectrum(cls, graph: np.ndarray) -> np.ndarray:
    """Sorted eigenvalues of the Laplacian L = D - A."""
    degrees = graph.sum(axis=1)
    L = np.diag(degrees) - graph
    return np.sort(np.linalg.eigvalsh(L))

get_edge_logit

get_edge_logit(sum_degrees: float) -> int

Bernoulli draw with probability = logistic(sum_degrees). Kept for backward compatibility; the add_remove_edge hot-path inlines a faster equivalent.

Source code in src/logit_graph/graph.py

def get_edge_logit(self, sum_degrees: float) -> int:
    """Bernoulli draw with probability = logistic(sum_degrees). Kept for backward
    compatibility; the add_remove_edge hot-path inlines a faster equivalent."""
    p = expit(sum_degrees)
    return int(self._rng.random() < p)

add_remove_edge

add_remove_edge() -> None

Propose a random node pair and set / clear their edge.

Source code in src/logit_graph/graph.py

def add_remove_edge(self) -> None:
    """Propose a random node pair and set / clear their edge."""
    i = int(self._rng.integers(0, self.n))
    j = int(self._rng.integers(0, self.n - 1))
    if j >= i:
        j += 1

    if self.layer2:
        feat = self._layer2_feature(i, j)
    else:
        sum_i = self._get_sum_degrees_fast(i)
        sum_j = self._get_sum_degrees_fast(j)
        if self.feature_mode == "paper_raw":
            feat = sum_i + sum_j
        elif self.feature_mode == "bounded":
            feat = np.log1p(sum_i) + np.log1p(sum_j)
        else:
            feat = self._layer2_feature(i, j)

    logit = self.sigma + self.alpha * self.beta * feat
    p = expit(logit)
    new_val = float(self._rng.random() < p)
    old_val = self._has_edge(i, j)
    self._sync_edge(i, j, new_val, old_val)

add_remove_edge_adj_only

add_remove_edge_adj_only() -> None

Neighbor-list path without dense-matrix sync (fast_mode partial).

Source code in src/logit_graph/graph.py

def add_remove_edge_adj_only(self) -> None:
    """Neighbor-list path without dense-matrix sync (fast_mode partial)."""
    prev = self._adj_only
    self._adj_only = True
    try:
        self.add_remove_edge()
    finally:
        self._adj_only = prev

materialize_adjacency

materialize_adjacency() -> None

Rebuild dense adjacency from neighbor lists.

Source code in src/logit_graph/graph.py

def materialize_adjacency(self) -> None:
    """Rebuild dense adjacency from neighbor lists."""
    n = self.n
    adj = np.zeros((n, n), dtype=float)
    for i in range(n):
        for j in self._nbrs[i]:
            adj[i, j] = 1.0
    self.graph = adj
    self._adj_only = False

populate_edges_baseline

populate_edges_baseline(
    warm_up: int,
    max_iterations: int,
    patience: int,
    check_interval: int = 50,
    edge_cv_tol: float = 0.02,
    spectrum_cv_tol: float = 0.02,
    fast_mode: bool = False,
) -> tuple[list[np.ndarray], np.ndarray]

Generate a graph without a ground-truth reference. Convergence uses the coefficient of variation (std/mean) of edge counts and spectrum norms over the last patience checks; fast_mode=True skips checks and runs max_iterations steps.

Source code in src/logit_graph/graph.py

def populate_edges_baseline(
    self,
    warm_up: int,
    max_iterations: int,
    patience: int,
    check_interval: int = 50,
    edge_cv_tol: float = 0.02,
    spectrum_cv_tol: float = 0.02,
    fast_mode: bool = False,
) -> tuple[list[np.ndarray], np.ndarray]:
    """Generate a graph without a ground-truth reference. Convergence uses the
    coefficient of variation (std/mean) of edge counts and spectrum norms over the last
    ``patience`` checks; ``fast_mode=True`` skips checks and runs ``max_iterations`` steps."""
    if fast_mode:
        fg = FastGibbsGraph(
            self.n,
            self.d,
            self.sigma,
            er_p=self.er_p,
            rng=self._rng,
            feature_mode=self.feature_mode,
            alpha=self.alpha,
            beta=self.beta,
            adj=self.graph,
        )
        fg.run_steps(max_iterations, self._rng)
        self._csr_rows = fg.to_rows_list()
        # Sweeps consume CSR rows directly; skip O(n²) adjacency + O(n³) spectrum.
        return [], np.empty(0, dtype=np.float64)

    graphs = deque(maxlen=max(patience + 10, 200))
    graphs.append(self.graph.copy())

    edge_history: deque[int] = deque(maxlen=patience)
    spectrum_norm_history: deque[float] = deque(maxlen=patience)

    pbar = tqdm(total=max_iterations, desc="Generating graph",
                leave=False, disable=False)

    for i in range(max_iterations):
        self.add_remove_edge()

        if i % check_interval == 0:
            graphs.append(self.graph.copy())
            edge_history.append(self._edge_count)

            spec_norm = float(np.linalg.norm(
                self.calculate_spectrum(self.graph)))
            spectrum_norm_history.append(spec_norm)

            if i >= warm_up and len(edge_history) >= patience:
                edges_arr = np.array(edge_history)
                spec_arr = np.array(spectrum_norm_history)

                mean_e = np.mean(edges_arr)
                cv_edges = (np.std(edges_arr) / mean_e) if mean_e > 0 else 0.0

                mean_s = np.mean(spec_arr)
                cv_spectrum = (np.std(spec_arr) / mean_s) if mean_s > 0 else 0.0

                stop_condition = (cv_edges < edge_cv_tol
                                  and cv_spectrum < spectrum_cv_tol)

                pbar.set_postfix({
                    'edges': self._edge_count,
                    'cv_e': f'{cv_edges:.4f}',
                    'cv_s': f'{cv_spectrum:.4f}',
                    'converged': stop_condition,
                })

                if stop_condition:
                    pbar.update(max_iterations - pbar.n)
                    break

        pbar.update(1)

    pbar.close()
    spectra = self.calculate_spectrum(self.graph)
    return list(graphs), spectra

populate_edges_spectrum

populate_edges_spectrum(
    warm_up: int,
    max_iterations: int,
    patience: int,
    real_graph: ndarray,
    edge_delta: Optional[float] = None,
    check_interval: int = 50,
    verbose: bool = True,
) -> tuple[list[np.ndarray], np.ndarray, list[float], int]

Legacy spectrum-only convergence (no GIC gate).

Source code in src/logit_graph/graph.py

def populate_edges_spectrum(
    self,
    warm_up: int,
    max_iterations: int,
    patience: int,
    real_graph: np.ndarray,
    edge_delta: Optional[float] = None,
    check_interval: int = 50,
    verbose: bool = True,
) -> tuple[list[np.ndarray], np.ndarray, list[float], int]:
    """Legacy spectrum-only convergence (no GIC gate)."""
    best_iteration = 0

    spectrum_diffs: list[float] = []
    real_spectrum = self.calculate_spectrum(real_graph)
    real_edges = int(np.triu(real_graph).sum())
    no_improvement_checks = 0
    best_spectrum_diff = float('inf')

    graphs = deque(maxlen=max(2 * patience + 100, 500))
    graphs.append(self.graph.copy())
    best_graph = self.graph.copy()

    for i in range(max_iterations):
        if edge_delta is not None:
            if self._edge_count > real_edges + edge_delta:
                if verbose:
                    print('Too many edges. Stopping.')
                break

        self.add_remove_edge()

        if i % check_interval == 0:
            graphs.append(self.graph.copy())
            current_spectrum = self.calculate_spectrum(self.graph)
            spectrum_diff = np.linalg.norm(current_spectrum - real_spectrum)
            spectrum_diffs.append(spectrum_diff)

            if verbose and i % 1000 == 0:
                print(f'\t Iteration {i}: spectrum diff = {spectrum_diff:.4f}')

            if spectrum_diff < best_spectrum_diff:
                best_spectrum_diff = spectrum_diff
                best_graph = self.graph.copy()
                best_iteration = i
                no_improvement_checks = 0
            elif i >= warm_up:
                no_improvement_checks += 1

            if i >= warm_up and no_improvement_checks >= patience:
                break

    if verbose:
        print(f'\t Best iteration: {best_iteration}')
        print(f'\t Best spectrum difference: {best_spectrum_diff:.4f}')
        print(f'\t Edges (current): {self._edge_count}, '
              f'Edges (real): {real_edges}')

    self.graph = best_graph
    self._init_cache()
    spectra = self.calculate_spectrum(self.graph)

    return list(graphs), spectra, spectrum_diffs, best_iteration

populate_edges_spectrum_min_gic

populate_edges_spectrum_min_gic(
    max_iterations: int,
    patience: int,
    real_graph: Union[ndarray, Graph],
    min_gic_threshold: float,
    gic_dist_type: str = "KL",
    edge_delta: Optional[float] = None,
    check_interval: int = 50,
    verbose: bool = True,
    er_p: float = 0.05,
) -> tuple[
    list[np.ndarray],
    np.ndarray,
    list[float],
    int,
    np.ndarray,
    list[float],
]

Populate edges targeting a real graph via a two-phase criterion: Phase 1 (GIC gate) iterates until GIC vs the real graph drops below min_gic_threshold; Phase 2 tracks the best Laplacian-spectrum distance and stops after patience checks w/o gain.

Source code in src/logit_graph/graph.py

def populate_edges_spectrum_min_gic(
    self,
    max_iterations: int,
    patience: int,
    real_graph: Union[np.ndarray, nx.Graph],
    min_gic_threshold: float,
    gic_dist_type: str = 'KL',
    edge_delta: Optional[float] = None,
    check_interval: int = 50,
    verbose: bool = True,
    er_p: float = 0.05,
) -> tuple[list[np.ndarray], np.ndarray, list[float], int, np.ndarray, list[float]]:
    """Populate edges targeting a real graph via a two-phase criterion: Phase 1 (GIC
    gate) iterates until GIC vs the real graph drops below ``min_gic_threshold``; Phase 2
    tracks the best Laplacian-spectrum distance and stops after ``patience`` checks w/o gain."""
    best_iteration = 0

    # GIC state
    gic_threshold_reached = False
    current_gic = float('inf')
    gic_values: list[float] = []

    # Spectrum state
    spectrum_diffs: list[float] = []
    if isinstance(real_graph, nx.Graph):
        real_graph_np = nx.to_numpy_array(real_graph)
    else:
        real_graph_np = real_graph

    real_spectrum = self.calculate_spectrum(real_graph_np)
    real_edges = int(np.triu(real_graph_np).sum())
    no_improvement_checks = 0
    best_spectrum_diff = float('inf')

    # Graph bookkeeping
    graphs = deque(maxlen=max(2 * patience + 100, 500))
    graphs.append(self.graph.copy())
    best_graph = self.graph.copy()

    # Pre-build NetworkX reference once for GIC calls
    real_nx_graph = nx.from_numpy_array(real_graph_np)

    # Progress bar
    pbar = None
    if verbose:
        pbar = tqdm(
            total=max_iterations,
            desc="Optimizing Graph", leave=True,
            bar_format=('{l_bar}{bar}| {n_fmt}/{total_fmt} '
                        '[{elapsed}<{remaining}, {rate_fmt}] {postfix}'))
        pbar.set_postfix({
            'GIC': f'{current_gic:.4f}',
            'Spec': f'{best_spectrum_diff:.4f}',
            'Pat': f'{no_improvement_checks}/{patience}',
            'Edges': f'{self._edge_count}/{real_edges}'
        })

    stop_reason = 'unknown'

    for i in range(max_iterations):
        # --- edge_delta guard (O(1) via cached count) ---
        if edge_delta is not None:
            if self._edge_count > real_edges + edge_delta:
                stop_reason = f'edge count exceeded delta ({edge_delta})'
                break

        # --- main step ---
        self.add_remove_edge()

        # --- periodic expensive check ---
        if i % check_interval == 0:
            graphs.append(self.graph.copy())

            current_spectrum = self.calculate_spectrum(self.graph)
            spectrum_diff = np.linalg.norm(current_spectrum - real_spectrum)
            spectrum_diffs.append(spectrum_diff)

            if not gic_threshold_reached:
                try:
                    current_nx = nx.from_numpy_array(self.graph)
                    gic_calc = gic.GraphInformationCriterion(
                        real_nx_graph, model='LG',
                        log_graph=current_nx, dist=gic_dist_type)
                    current_gic = gic_calc.calculate_spectral_distance()
                except Exception:
                    current_gic = float('inf')

                if current_gic <= min_gic_threshold:
                    if verbose and pbar:
                        pbar.write(
                            f'GIC threshold {min_gic_threshold} reached '
                            f'at iteration {i:,} (GIC: {current_gic:.4f}). '
                            f'Starting spectrum patience ({patience} checks).')
                    gic_threshold_reached = True
                    no_improvement_checks = 0

            gic_values.append(current_gic)

            if spectrum_diff < best_spectrum_diff:
                best_spectrum_diff = spectrum_diff
                best_graph = self.graph.copy()
                best_iteration = i
                if gic_threshold_reached:
                    no_improvement_checks = 0
            elif gic_threshold_reached:
                no_improvement_checks += 1

            if verbose and pbar:
                pbar.set_postfix({
                    'GIC': f'{current_gic:.4f}',
                    'Spec': f'{best_spectrum_diff:.4f}',
                    'Pat': f'{no_improvement_checks}/{patience}',
                    'Edges': f'{self._edge_count}/{real_edges}'
                })

            if gic_threshold_reached and no_improvement_checks >= patience:
                stop_reason = (
                    f'no spectral improvement for {patience} checks '
                    f'after GIC threshold was met')
                break

        if verbose and pbar:
            pbar.update(1)
    else:
        stop_reason = f'max iterations ({max_iterations:,}) reached'

    # --- summary (before closing pbar) ---
    if verbose and pbar:
        pbar.write(f'\nStopping: {stop_reason}')
        pbar.write(f'  Best iteration: {best_iteration:,}')
        pbar.write(f'  Best spectrum diff: {best_spectrum_diff:.4f}')
        pbar.write(f'  Edges in best graph: '
                   f'{int(np.triu(best_graph).sum())} '
                   f'(real: {real_edges})')
        pbar.close()

    spectra = self.calculate_spectrum(best_graph)
    return (list(graphs), spectra, spectrum_diffs,
            best_iteration, best_graph, gic_values)

LogitRegEstimator

LogitRegEstimator(
    graph: Union[ndarray, Graph],
    d: int,
    verbose: bool = False,
    layer2: bool = True,
    feature_mode: FeatureMode = "incremental",
)

Source code in src/logit_graph/logit_estimator.py

def __init__(
    self,
    graph: Union[np.ndarray, nx.Graph],
    d: int,
    verbose: bool = False,
    layer2: bool = True,
    feature_mode: FeatureMode = "incremental",
) -> None:
    # NOTE: feature_mode defaults to "incremental" to match the generator and the
    # paper's offset-logit form; the legacy "bounded" mode adds a nonzero pair term
    # even at d=0, biasing sigma_hat away from logit(empirical_density).
    self.graph = graph  # The observed adjacency matrix
    if isinstance(graph, np.ndarray):
        self.n = graph.shape[0]  # Number of nodes in the graph
        self.adj = graph
    elif isinstance(graph, nx.Graph):
        self.n = graph.number_of_nodes()
        self.adj = nx.to_numpy_array(graph)
    else:
        raise ValueError("Unsupported graph type. Please provide a NumPy array or NetworkX graph.")
    self.d = d # number of degrees to search
    self.verbose = verbose
    self.layer2 = layer2
    self.feature_mode = feature_mode

compute_aic

compute_aic(
    d_est: Optional[int] = None,
    feature_mode: Optional[FeatureMode] = None,
    layer2: Optional[bool] = None,
    extra_penalty: float = 0.0,
) -> dict[str, float]

Paper AIC: -2*ll + 2 (one parameter sigma).

Source code in src/logit_graph/logit_estimator.py

def compute_aic(
    self,
    d_est: Optional[int] = None,
    feature_mode: Optional[FeatureMode] = None,
    layer2: Optional[bool] = None,
    extra_penalty: float = 0.0,
) -> dict[str, float]:
    """Paper AIC: -2*ll + 2 (one parameter sigma)."""
    d_use = self.d if d_est is None else d_est
    use_mode = self.feature_mode if feature_mode is None else feature_mode
    use_layer2 = self.layer2 if layer2 is None else layer2

    offsets, labels_arr = build_pair_dataset(
        self.adj,
        d=d_use,
        mode=use_mode,
        layer2=use_layer2,
    )
    from .offset_logit import aic_from_offset_fit, fit_offset_logit_fast

    sigma_hat, ll = fit_offset_logit_fast(offsets, labels_arr)
    fit = aic_from_offset_fit(sigma_hat, ll, extra_penalty=extra_penalty)
    return {
        "aic": fit["aic"],
        "ll": fit["ll"],
        "k": fit["k"],
        "sigma_hat": fit["sigma_hat"],
        "d_est": float(d_use),
        "n_obs": float(len(labels_arr)),
    }

select_d

select_d(
    d_candidates: Optional[list[int]] = None,
    feature_mode: Optional[FeatureMode] = None,
    extra_penalty_per_d: float = 0.0,
) -> tuple[int, dict[int, dict[str, float]]]

Return argmin_d AIC(d) and per-d stats.

Source code in src/logit_graph/logit_estimator.py

def select_d(
    self,
    d_candidates: Optional[list[int]] = None,
    feature_mode: Optional[FeatureMode] = None,
    extra_penalty_per_d: float = 0.0,
) -> tuple[int, dict[int, dict[str, float]]]:
    """Return argmin_d AIC(d) and per-d stats."""
    if d_candidates is None:
        d_candidates = [0, 1, 2, 3]
    use_mode = self.feature_mode if feature_mode is None else feature_mode
    stats: dict[int, dict[str, float]] = {}
    for d_c in d_candidates:
        stats[d_c] = self.compute_aic(
            d_est=d_c,
            feature_mode=use_mode,
            extra_penalty=extra_penalty_per_d * d_c,
        )
    best = min(stats, key=lambda d: stats[d]["aic"])
    return best, stats

estimate_parameters

estimate_parameters(
    l1_wt: float = 1,
    alpha: float = 0,
    features: Optional[ndarray] = None,
    labels: Optional[list[int]] = None,
    fix_beta: bool = True,
) -> tuple[Any, np.ndarray, np.ndarray]

Fit the logistic regression on the extracted features/labels. fix_beta=True (default) fixes the degree-sum coefficient at 1 (offset) and estimates only sigma (paper form); fix_beta=False frees both sigma and beta, using l1_wt/alpha (legacy).

Source code in src/logit_graph/logit_estimator.py

def estimate_parameters(
    self,
    l1_wt: float = 1,
    alpha: float = 0,
    features: Optional[np.ndarray] = None,
    labels: Optional[list[int]] = None,
    fix_beta: bool = True,
) -> tuple[Any, np.ndarray, np.ndarray]:
    """Fit the logistic regression on the extracted features/labels. fix_beta=True
    (default) fixes the degree-sum coefficient at 1 (offset) and estimates only sigma
    (paper form); fix_beta=False frees both sigma and beta, using l1_wt/alpha (legacy)."""
    if features is None or labels is None:
        features, labels = self.get_features_labels()

    if self.verbose:
        print("\nStarting parameter estimation...")
        print(f"fix_beta={fix_beta}")

    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", message="overflow encountered in exp")
        warnings.filterwarnings("ignore", message="divide by zero encountered in log")
        warnings.filterwarnings("ignore", category=RuntimeWarning, module="statsmodels")

        if fix_beta:
            intercept_col = features[:, 0:1]
            offset_col = features[:, 1]
            model = sm.Logit(labels, intercept_col, offset=offset_col)
            try:
                result = model.fit(method='bfgs', disp=self.verbose)
            except (np.linalg.LinAlgError, Exception):
                result = model.fit_regularized(
                    method='l1', alpha=1e-4, disp=self.verbose,
                )
        else:
            model = sm.Logit(labels, features)
            if l1_wt in [0, 1]:
                result = model.fit_regularized(
                    method='l1' if l1_wt == 1 else None,
                    alpha=alpha, disp=self.verbose,
                )
            else:
                result = model.fit_regularized(
                    method='elastic_net', alpha=alpha,
                    L1_wt=l1_wt, disp=self.verbose,
                )

    params = result.params

    if self.verbose:
        param_vals = params.values if hasattr(params, 'values') else params
        print(f"Estimated parameters: {param_vals}")
        mle_ret = getattr(result, 'mle_retvals', None)
        if mle_ret:
            print(f"Converged: {mle_ret.get('converged', 'N/A')}")

    return result, params, features

Features

build_pair_dataset

build_pair_dataset(
    graph: Union[ndarray, Graph],
    d: int,
    mode: FeatureMode = "bounded",
    layer2: bool = True,
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
    max_pairs: Optional[int] = None,
    seed: Optional[int] = None,
) -> tuple[np.ndarray, np.ndarray]

Build offset vector and binary labels for all upper-triangle pairs. Returns (offsets, labels): offsets are (m,) feature values (the fixed beta=1 offset), labels are (m,) 0/1 edge indicators.

Source code in src/logit_graph/lg_features.py

def build_pair_dataset(
    graph: Union[np.ndarray, nx.Graph],
    d: int,
    mode: FeatureMode = "bounded",
    layer2: bool = True,
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
    max_pairs: Optional[int] = None,
    seed: Optional[int] = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Build offset vector and binary labels for all upper-triangle pairs. Returns
    ``(offsets, labels)``: offsets are (m,) feature values (the fixed beta=1 offset),
    labels are (m,) 0/1 edge indicators."""
    adj = _adj_from_input(graph)
    n = adj.shape[0]

    if layer2 and max_pairs is None:
        from .lg_features_fast import build_pair_dataset_fast

        return build_pair_dataset_fast(adj, d, mode=mode, alpha_gwesp=alpha_gwesp)

    rng = np.random.default_rng(seed)

    pairs: list[tuple[int, int]] = []
    for i in range(n):
        for j in range(i + 1, n):
            pairs.append((i, j))

    if max_pairs is not None and len(pairs) > max_pairs:
        idx = rng.choice(len(pairs), size=max_pairs, replace=False)
        pairs = [pairs[int(k)] for k in idx]

    use_vertex_cache = layer2 and mode in ("paper_raw", "bounded")
    vertex_sums: Optional[np.ndarray] = None
    if use_vertex_cache:
        vertex_sums = precompute_vertex_sums(adj, d)

    offsets = np.empty(len(pairs), dtype=float)
    labels = np.empty(len(pairs), dtype=int)

    for k, (i, j) in enumerate(pairs):
        labels[k] = 1 if adj[i, j] > 0 else 0
        if layer2:
            if use_vertex_cache and vertex_sums is not None:
                adj_l2 = adj
                si = float(vertex_sums[i])
                sj = float(vertex_sums[j])
                if adj[i, j] > 0:
                    si -= 1.0
                    sj -= 1.0
                if mode == "paper_raw":
                    offsets[k] = si + sj
                else:
                    offsets[k] = math.log(1.0 + si) + math.log(1.0 + sj)
            else:
                offsets[k] = pair_feature_layer2(
                    adj, i, j, d, mode=mode, alpha_gwesp=alpha_gwesp,
                )
        else:
            offsets[k] = pair_feature(
                adj, i, j, d, mode=mode, alpha_gwesp=alpha_gwesp,
            )

    return offsets, labels

pair_feature

pair_feature(
    graph: ndarray,
    i: int,
    j: int,
    d: int,
    mode: FeatureMode = "bounded",
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
) -> float

Pair feature x_ij used as offset (beta fixed at 1 in paper estimator).

Source code in src/logit_graph/lg_features.py

def pair_feature(
    graph: np.ndarray,
    i: int,
    j: int,
    d: int,
    mode: FeatureMode = "bounded",
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
) -> float:
    """Pair feature x_ij used as offset (beta fixed at 1 in paper estimator)."""
    if mode == "paper_raw":
        return sum_degree(graph, i, d) + sum_degree(graph, j, d)
    if mode == "bounded":
        si = sum_degree(graph, i, d)
        sj = sum_degree(graph, j, d)
        return math.log(1.0 + si) + math.log(1.0 + sj)
    if mode == "incremental":
        return incremental_h(graph, i, j, d, alpha_gwesp=alpha_gwesp)
    if mode == "common_dhop":
        return math.log(1.0 + common_dhop_count(graph, i, j, d))
    raise ValueError(f"Unknown feature mode: {mode!r}")

pair_feature_layer2

pair_feature_layer2(
    graph: Union[ndarray, Graph],
    i: int,
    j: int,
    d: int,
    mode: FeatureMode = "bounded",
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
) -> float

Leave-one-out (full conditional) feature: computed on the graph with edge (i,j) removed.

Source code in src/logit_graph/lg_features.py

def pair_feature_layer2(
    graph: Union[np.ndarray, nx.Graph],
    i: int,
    j: int,
    d: int,
    mode: FeatureMode = "bounded",
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
) -> float:
    """Leave-one-out (full conditional) feature: computed on the graph with edge (i,j) removed."""
    adj = _adj_from_input(graph).copy()
    had = adj[i, j] > 0
    if had:
        adj[i, j] = adj[j, i] = 0.0
    return pair_feature(adj, i, j, d, mode=mode, alpha_gwesp=alpha_gwesp)

incremental_h

incremental_h(
    graph: ndarray,
    i: int,
    j: int,
    d: int,
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
) -> float

Source code in src/logit_graph/lg_features.py

def incremental_h(
    graph: np.ndarray,
    i: int,
    j: int,
    d: int,
    alpha_gwesp: float = ALPHA_GWESP_DEFAULT,
) -> float:
    if d == 0:
        return 0.0
    if d == 1:
        c = common_dhop_count(graph, i, j, 1)
        if c <= 0:
            return 0.0
        return alpha_gwesp * (1.0 - (1.0 - 1.0 / alpha_gwesp) ** c)
    delta = common_dhop_count(graph, i, j, d) - common_dhop_count(graph, i, j, d - 1)
    return math.log(1.0 + max(0, delta))

recommended_iterations

recommended_iterations(
    n: int, cap: Optional[int] = None
) -> int

Gibbs steps so each edge is resampled O(1) times.

Source code in src/logit_graph/lg_features.py

def recommended_iterations(n: int, cap: Optional[int] = None) -> int:
    """Gibbs steps so each edge is resampled O(1) times."""
    base = max(20_000, int(5 * n * (n - 1)))
    if cap is not None:
        return min(base, cap)
    return base

LG: iterative generation & estimation

The LG iterative generation algorithm and its estimator: logit Pr[A_ij(t)=1 | G(t-1)] = σ + α·D(t-1) (Eq. 3.5), with the degree term D = g(S_i)+g(S_j) read from the predetermined previous graph so the pooled dyad design is an ordinary logistic regression and the MLE of (σ, α) is consistent.

grow_graph

grow_graph(
    n: int,
    d: int,
    sigma: float,
    alpha: float,
    *,
    n_steps: int,
    degree_mode: FeatureMode = DEGREE_MODE,
    seed: int | None = None,
    p0: float = 0.02,
    record_design: bool = True,
    store_snapshots: bool = True,
    allow_removal: bool = True,
    until_convergence: bool = False,
    esd_tol: float = 0.01,
    patience: int = 3,
    esd_nbins: int = 50,
    step_callback=None
) -> GrowthResult

Grow a temporal Logit-Graph from a sparse ER seed (degree-only). allow_removal (default) resamples every dyad from lagged expit(sigma+alpha*D) — edges form and dissolve (ergodic chain); allow_removal=False grows add-only to saturation.

Source code in src/logit_graph/temporal.py

def grow_graph(
    n: int,
    d: int,
    sigma: float,
    alpha: float,
    *,
    n_steps: int,
    degree_mode: FeatureMode = DEGREE_MODE,
    seed: int | None = None,
    p0: float = 0.02,
    record_design: bool = True,
    store_snapshots: bool = True,
    allow_removal: bool = True,
    until_convergence: bool = False,
    esd_tol: float = 1e-2,
    patience: int = 3,
    esd_nbins: int = 50,
    step_callback=None,
) -> GrowthResult:
    """Grow a temporal Logit-Graph from a sparse ER seed (degree-only). allow_removal
    (default) resamples every dyad from lagged expit(sigma+alpha*D) — edges form and
    dissolve (ergodic chain); allow_removal=False grows add-only to saturation."""
    rng = np.random.default_rng(seed)
    rows, cols = np.triu_indices(n, k=1)

    adj = np.zeros((n, n), dtype=np.float64)
    seed_mask = rng.random(rows.shape[0]) < p0
    adj[rows[seed_mask], cols[seed_mask]] = 1.0
    adj[cols[seed_mask], rows[seed_mask]] = 1.0

    snapshots = [adj.copy()] if store_snapshots else []
    Xs: list[np.ndarray] = []
    ys: list[np.ndarray] = []

    esd_kl_trace: list[float] = []
    below = 0
    converged = False
    n_steps_run = n_steps

    for step in range(n_steps):
        prev_adj = adj.copy() if until_convergence else None
        D, labels = _degree_feature(adj, d, degree_mode)
        p = expit(sigma + alpha * D)
        draw = rng.random(p.shape[0]) < p
        if allow_removal:
            # Resample every dyad from the lagged probability: edges may form OR
            # dissolve. Design is over all dyads, outcome = the new state.
            if record_design:
                Xs.append(D.reshape(-1, 1))
                ys.append(draw.astype(np.int8))
            adj[:] = 0.0
            ki, kj = rows[draw], cols[draw]
            adj[ki, kj] = 1.0
            adj[kj, ki] = 1.0
        else:
            at_risk = labels == 0
            form = at_risk & draw
            if record_design:
                Xs.append(D[at_risk].reshape(-1, 1))
                ys.append(form[at_risk].astype(np.int8))
            fi, fj = rows[form], cols[form]
            adj[fi, fj] = 1.0
            adj[fj, fi] = 1.0
        if store_snapshots:
            snapshots.append(adj.copy())

        if until_convergence:
            kl = _adjacency_esd_kl(adj, prev_adj, esd_nbins)
            esd_kl_trace.append(kl)
            below = below + 1 if kl < esd_tol else 0
            if below >= patience:
                converged = True
                n_steps_run = step + 1
                break

        # Generic per-step observer: receives (step, live adj) — copy if retaining.
        # Returning True stops growth early (e.g. GIC early-stopping by a caller).
        if step_callback is not None and step_callback(step, adj):
            n_steps_run = step + 1
            break

    X = np.vstack(Xs) if Xs else np.empty((0, 1), dtype=np.float64)
    y = np.concatenate(ys) if ys else np.empty(0, dtype=np.int8)
    params = dict(n=n, d=d, sigma=sigma, alpha=alpha, n_steps=n_steps,
                  degree_mode=degree_mode, p0=p0, allow_removal=allow_removal,
                  n_steps_run=n_steps_run)
    if until_convergence:
        params.update(until_convergence=True, esd_tol=esd_tol, patience=patience,
                      n_steps_run=n_steps_run, converged=converged,
                      esd_kl_trace=esd_kl_trace)
    return GrowthResult(adj=adj, X=X, y=y, snapshots=snapshots, params=params)

GrowthResult dataclass

GrowthResult(
    adj: ndarray,
    X: ndarray,
    y: ndarray,
    snapshots: list,
    params: dict,
)

Output of :func:grow_graph: final adjacency adj (n x n, 0/1, symmetric), pooled design X (m,1)=[D] / outcomes y (1=at-risk non-edge formed), snapshots G(0..n_steps) (empty if not stored), and the generative params.

growth_design_from_snapshots

growth_design_from_snapshots(
    snapshots: list,
    d: int,
    *,
    degree_mode: FeatureMode = DEGREE_MODE,
    allow_removal: bool = True
) -> tuple[np.ndarray, np.ndarray]

Build the at-risk logistic-regression design from observed snapshots: predictor D from G(t-1), outcome = formed in G(t). With allow_removal the set is all dyads and the outcome is each dyad's full new state — matching grow_graph (whose design it reproduces).

Source code in src/logit_graph/temporal.py

def growth_design_from_snapshots(
    snapshots: list,
    d: int,
    *,
    degree_mode: FeatureMode = DEGREE_MODE,
    allow_removal: bool = True,
) -> tuple[np.ndarray, np.ndarray]:
    """Build the at-risk logistic-regression design from observed snapshots: predictor D
    from G(t-1), outcome = formed in G(t). With allow_removal the set is all dyads and the
    outcome is each dyad's full new state — matching grow_graph (whose design it reproduces)."""
    Xs: list[np.ndarray] = []
    ys: list[np.ndarray] = []
    for t in range(len(snapshots) - 1):
        prev = np.asarray(snapshots[t])
        nxt = np.asarray(snapshots[t + 1])
        n = prev.shape[0]
        rows, cols = np.triu_indices(n, k=1)
        D, labels_prev = _degree_feature(prev, d, degree_mode)
        state = (nxt[rows, cols] > 0)
        if allow_removal:
            Xs.append(D.reshape(-1, 1))
            ys.append(state.astype(np.int8))
        else:
            at_risk = labels_prev == 0
            Xs.append(D[at_risk].reshape(-1, 1))
            ys.append(state[at_risk].astype(np.int8))
    X = np.vstack(Xs) if Xs else np.empty((0, 1), dtype=np.float64)
    y = np.concatenate(ys) if ys else np.empty(0, dtype=np.int8)
    return X, y

fit_growth_params

fit_growth_params(X: ndarray, labels: ndarray) -> dict

Fit logit(P[form]) = sigma + alpha*D by ordinary logistic regression — the exact, consistent MLE for the growth model (dyad-independent given the past). Returns sigma, alpha, their SEs, the log-likelihood, AIC (k=2), and n_params.

Source code in src/logit_graph/temporal.py

def fit_growth_params(X: np.ndarray, labels: np.ndarray) -> dict:
    """Fit logit(P[form]) = sigma + alpha*D by ordinary logistic regression — the exact,
    consistent MLE for the growth model (dyad-independent given the past). Returns sigma,
    alpha, their SEs, the log-likelihood, AIC (k=2), and n_params."""
    y = np.asarray(labels, dtype=int)
    Xc = sm.add_constant(np.asarray(X, dtype=np.float64), has_constant="add")
    res = None
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        for method in ("newton", "bfgs", "lbfgs"):
            try:
                r = sm.Logit(y, Xc).fit(method=method, disp=0, maxiter=200)
                if np.isfinite(r.llf):
                    res = r
                    break
            except Exception:
                continue
        if res is None:
            res = sm.Logit(y, Xc).fit_regularized(method="l1", alpha=1e-4, disp=0)

    params = np.asarray(res.params, dtype=np.float64)
    try:
        bse = np.asarray(res.bse, dtype=np.float64)
    except Exception:
        bse = np.full(2, np.nan)
    ll = float(res.llf)
    return {
        "sigma": float(params[0]),
        "alpha": float(params[1]),
        "se_sigma": float(bse[0]),
        "se_alpha": float(bse[1]),
        "ll": ll,
        "aic": float(aic_from_offset_fit(params[0], ll, k=2.0)["aic"]),
        "n_params": 2,
    }

fit_growth_from_result

fit_growth_from_result(result: GrowthResult) -> dict

Convenience: fit on the design recorded by :func:grow_graph.

Source code in src/logit_graph/temporal.py

def fit_growth_from_result(result: GrowthResult) -> dict:
    """Convenience: fit on the design recorded by :func:`grow_graph`."""
    return fit_growth_params(result.X, result.y)

Unified LG: general pairwise features

The §3.7 extension logit p_ij = σ + Σₖ θₖ·φₖ(i,j): the degree term plus arbitrary appended pairwise features. grow_graph_multi / fit_multi_params mirror the degree-only API but carry the extra features; community_feature (block membership) and latent_feature (latent proximity) build the canonical covariates of the thesis's unified five-feature LG model from an observed graph.

grow_graph_multi

grow_graph_multi(
    n: int,
    d: int,
    sigma: float,
    alpha: float,
    features: ndarray,
    coefs,
    *,
    n_steps: int,
    feature_names: list | None = None,
    degree_mode: FeatureMode = DEGREE_MODE,
    allow_removal: bool = True,
    seed: int | None = None,
    p0: float = 0.02,
    record_design: bool = True,
    store_snapshots: bool = True
) -> MultiGrowthResult

Generate from the unified LG (general pairwise features) by growth from a sparse ER seed: features (n(n-1)/2, k) are fixed pairwise covariates phi_k with coefficients coefs, and each step draws dyads from expit(sigma + alphaD + features @ coefs), D from G(t-1) (allow_removal=ergodic).

Source code in src/logit_graph/temporal_multi.py

def grow_graph_multi(
    n: int,
    d: int,
    sigma: float,
    alpha: float,
    features: np.ndarray,
    coefs,
    *,
    n_steps: int,
    feature_names: list | None = None,
    degree_mode: FeatureMode = DEGREE_MODE,
    allow_removal: bool = True,
    seed: int | None = None,
    p0: float = 0.02,
    record_design: bool = True,
    store_snapshots: bool = True,
) -> MultiGrowthResult:
    """Generate from the unified LG (general pairwise features) by growth from a sparse ER seed:
    ``features`` (n*(n-1)/2, k) are fixed pairwise covariates phi_k with coefficients ``coefs``, and each
    step draws dyads from expit(sigma + alpha*D + features @ coefs), D from G(t-1) (allow_removal=ergodic)."""
    features = np.atleast_2d(np.asarray(features, dtype=np.float64))
    if features.shape[0] != n * (n - 1) // 2:
        features = features.T
    coefs = np.asarray(coefs, dtype=np.float64).ravel()
    k = features.shape[1]
    if coefs.shape[0] != k:
        raise ValueError(f"coefs has length {coefs.shape[0]} but features has {k} columns")
    if feature_names is None:
        feature_names = [f"f{j}" for j in range(k)]
    fixed_lo = features @ coefs

    rng = np.random.default_rng(seed)
    rows, cols = np.triu_indices(n, k=1)
    adj = np.zeros((n, n), dtype=np.float64)
    seed_mask = rng.random(rows.shape[0]) < p0
    adj[rows[seed_mask], cols[seed_mask]] = 1.0
    adj[cols[seed_mask], rows[seed_mask]] = 1.0

    snapshots = [adj.copy()] if store_snapshots else []
    Xs: list[np.ndarray] = []
    ys: list[np.ndarray] = []

    for _ in range(n_steps):
        D, labels = _degree_feature(adj, d, degree_mode)
        p = expit(sigma + alpha * D + fixed_lo)
        draw = rng.random(p.shape[0]) < p
        design = np.column_stack([D, features])
        if allow_removal:
            if record_design:
                Xs.append(design)
                ys.append(draw.astype(np.int8))
            adj[:] = 0.0
            adj[rows[draw], cols[draw]] = 1.0
            adj[cols[draw], rows[draw]] = 1.0
        else:
            at_risk = labels == 0
            form = at_risk & draw
            if record_design:
                Xs.append(design[at_risk])
                ys.append(form[at_risk].astype(np.int8))
            adj[rows[form], cols[form]] = 1.0
            adj[cols[form], rows[form]] = 1.0
        if store_snapshots:
            snapshots.append(adj.copy())

    X = np.vstack(Xs) if Xs else np.empty((0, 1 + k), dtype=np.float64)
    y = np.concatenate(ys) if ys else np.empty(0, dtype=np.int8)
    params = dict(n=n, d=d, sigma=sigma, alpha=alpha, coefs=coefs.tolist(),
                  feature_names=list(feature_names), n_steps=n_steps,
                  degree_mode=degree_mode, allow_removal=allow_removal, p0=p0)
    return MultiGrowthResult(adj=adj, X=X, y=y, feature_names=list(feature_names),
                             snapshots=snapshots, params=params)

MultiGrowthResult dataclass

MultiGrowthResult(
    adj: ndarray,
    X: ndarray,
    y: ndarray,
    feature_names: list,
    snapshots: list,
    params: dict,
)

Output of :func:grow_graph_multi: final adjacency adj (n x n, 0/1, symmetric), pooled design X (m, 1+k) = [D, F_1..F_k] / outcomes y, the feature_names for the extra covariates, the snapshots G(0..n_steps) (empty if not stored), and params.

multi_design_from_snapshots

multi_design_from_snapshots(
    snapshots: list,
    d: int,
    features: ndarray,
    *,
    degree_mode: FeatureMode = DEGREE_MODE,
    allow_removal: bool = True
) -> tuple[np.ndarray, np.ndarray]

Build the pooled [D, F_1..F_k] design from observed snapshots: D from G(t-1), the extra fixed features repeated each step, outcome from G(t). With allow_removal the design is over all dyads (else the at-risk non-edges). Reproduces grow_graph_multi's recorded design.

Source code in src/logit_graph/temporal_multi.py

def multi_design_from_snapshots(
    snapshots: list,
    d: int,
    features: np.ndarray,
    *,
    degree_mode: FeatureMode = DEGREE_MODE,
    allow_removal: bool = True,
) -> tuple[np.ndarray, np.ndarray]:
    """Build the pooled [D, F_1..F_k] design from observed snapshots: D from G(t-1), the extra
    fixed ``features`` repeated each step, outcome from G(t). With allow_removal the design is
    over all dyads (else the at-risk non-edges). Reproduces grow_graph_multi's recorded design."""
    features = np.atleast_2d(np.asarray(features, dtype=np.float64))
    Xs: list[np.ndarray] = []
    ys: list[np.ndarray] = []
    for t in range(len(snapshots) - 1):
        prev = np.asarray(snapshots[t])
        nxt = np.asarray(snapshots[t + 1])
        n = prev.shape[0]
        if features.shape[0] != n * (n - 1) // 2:
            features = features.T
        rows, cols = np.triu_indices(n, k=1)
        D, labels_prev = _degree_feature(prev, d, degree_mode)
        state = (nxt[rows, cols] > 0)
        design = np.column_stack([D, features])
        if allow_removal:
            Xs.append(design)
            ys.append(state.astype(np.int8))
        else:
            at_risk = labels_prev == 0
            Xs.append(design[at_risk])
            ys.append(state[at_risk].astype(np.int8))
    k = features.shape[1]
    X = np.vstack(Xs) if Xs else np.empty((0, 1 + k), dtype=np.float64)
    y = np.concatenate(ys) if ys else np.empty(0, dtype=np.int8)
    return X, y

fit_multi_params

fit_multi_params(
    X: ndarray,
    labels: ndarray,
    feature_names: list | None = None,
) -> dict

Fit logit(P) = sigma + alphaD + sum_k beta_kF_k by pooled logistic regression — the exact, consistent MLE for the multi-feature temporal model. Returns sigma, alpha, the per-feature coefficients and their SEs, the log-likelihood, AIC, and n_params (solver fallbacks for robustness).

Source code in src/logit_graph/temporal_multi.py

def fit_multi_params(X: np.ndarray, labels: np.ndarray, feature_names: list | None = None) -> dict:
    """Fit logit(P) = sigma + alpha*D + sum_k beta_k*F_k by pooled logistic regression — the exact,
    consistent MLE for the multi-feature temporal model. Returns sigma, alpha, the per-feature
    coefficients and their SEs, the log-likelihood, AIC, and n_params (solver fallbacks for robustness)."""
    y = np.asarray(labels, dtype=int)
    X = np.asarray(X, dtype=np.float64)
    k = X.shape[1] - 1
    if feature_names is None:
        feature_names = [f"f{j}" for j in range(k)]
    Xc = sm.add_constant(X, has_constant="add")
    res = None
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        for method in ("newton", "bfgs", "lbfgs"):
            try:
                r = sm.Logit(y, Xc).fit(method=method, disp=0, maxiter=300)
                if np.isfinite(r.llf):
                    res = r
                    break
            except Exception:
                continue
        if res is None:
            res = sm.Logit(y, Xc).fit_regularized(method="l1", alpha=1e-4, disp=0)

    b = np.asarray(res.params, dtype=np.float64)
    try:
        se = np.asarray(res.bse, dtype=np.float64)
    except Exception:
        se = np.full(b.shape[0], np.nan)
    ll = float(res.llf)
    n_params = 2 + k
    out = {
        "sigma": float(b[0]),
        "alpha": float(b[1]),
        "se_sigma": float(se[0]),
        "se_alpha": float(se[1]),
        "coefs": {feature_names[j]: float(b[2 + j]) for j in range(k)},
        "se_coefs": {feature_names[j]: float(se[2 + j]) for j in range(k)},
        "ll": ll,
        "aic": 2.0 * n_params - 2.0 * ll,
        "n_params": n_params,
    }
    return out

community_feature

community_feature(
    graph,
    *,
    resolution: float = 1.0,
    seed: int | None = None
) -> np.ndarray

Same-community indicator over upper-triangle dyads from a fixed Louvain partition of the graph at the given resolution (higher resolution -> finer communities). An exogenous node covariate, like SBM uses, so it stays identifiable. Returns a (n*(n-1)/2,) float array.

Source code in src/logit_graph/temporal_multi.py

def community_feature(graph, *, resolution: float = 1.0, seed: int | None = None) -> np.ndarray:
    """Same-community indicator over upper-triangle dyads from a fixed Louvain partition of the
    graph at the given resolution (higher resolution -> finer communities). An exogenous node
    covariate, like SBM uses, so it stays identifiable. Returns a (n*(n-1)/2,) float array."""
    import networkx as nx
    A = _as_adj(graph)
    n = A.shape[0]
    rows, cols = np.triu_indices(n, k=1)
    part = nx.community.louvain_communities(nx.from_numpy_array(A), seed=seed, resolution=resolution)
    blk = np.empty(n, dtype=int)
    for i, com in enumerate(part):
        for v in com:
            blk[v] = i
    return (blk[rows] == blk[cols]).astype(np.float64)

latent_feature

latent_feature(
    graph, k: int = 4, *, kind: str = "dot"
) -> np.ndarray

Latent proximity over upper-triangle dyads from a rank-k adjacency spectral embedding z (eigvecs scaled by sqrt|eigval|). kind="dot": z_i . z_j (RDPG inner product); kind="dist": -||z_i - z_j|| (Hoff latent-space distance). Standardized; returns a (n*(n-1)/2,) array.

Source code in src/logit_graph/temporal_multi.py

def latent_feature(graph, k: int = 4, *, kind: str = "dot") -> np.ndarray:
    """Latent proximity over upper-triangle dyads from a rank-k adjacency spectral embedding z
    (eigvecs scaled by sqrt|eigval|). kind="dot": z_i . z_j (RDPG inner product); kind="dist":
    -||z_i - z_j|| (Hoff latent-space distance). Standardized; returns a (n*(n-1)/2,) array."""
    A = _as_adj(graph)
    n = A.shape[0]
    rows, cols = np.triu_indices(n, k=1)
    w, U = np.linalg.eigh(A)
    idx = np.argsort(-np.abs(w))[:k]
    z = U[:, idx] * np.sqrt(np.abs(w[idx]))
    if kind == "dist":
        L = -np.sqrt(((z[rows] - z[cols]) ** 2).sum(1))
    else:
        L = (z[rows] * z[cols]).sum(1)
    return (L - L.mean()) / (L.std() + 1e-9)

Experiment configuration

AICSweepConfig dataclass

AICSweepConfig(
    d_true_values: list[int] = (lambda: [0, 1, 2, 3])(),
    d_est_values: list[int] = (lambda: [0, 1, 2, 3])(),
    n_sizes: list[int] = (lambda: [100])(),
    n_runs: int = 4,
    m_ensemble: int = 3,
    iter_cap: Optional[int] = 30000,
    target_density: float = 0.1,
    signal: float = 0.5,
    aic_penalty_per_d: float = 3.0,
    feature_mode_gen: str = "incremental",
    feature_mode_est: str = "incremental",
    seed_base: int = 1000,
    sigma_gen: float = -3.0,
    iter_cap_by_d: Optional[dict] = None,
    sigma_gen_per_n: Optional[dict] = None,
    m_ensemble_by_d: Optional[dict] = None,
    sigma_gen_per_n_d: Optional[dict] = None,
)

SigmaSweepConfig dataclass

SigmaSweepConfig(
    sigma_values: list[float] = (
        lambda: [-2.0, -4.0, -6.0, -8.0]
    )(),
    d_values: list[int] = (lambda: [0, 1, 2])(),
    n_values: list[int] = (lambda: [50, 100])(),
    n_reps: int = 4,
    iter_cap: Optional[int] = 80000,
    target_density: float = 0.1,
    signal: float = 0.5,
    feature_mode_gen: str = "incremental",
    feature_mode_est: str = "incremental",
    seed_base: int = 0,
    adaptive_stopping: bool = False,
    adaptive_check_interval: int = 20000,
    adaptive_patience: int = 3,
    adaptive_cv_tol: float = 0.02,
    adaptive_min_iter: int = 20000,
    iter_cap_by_d: Optional[dict] = None,
    n_values_by_sigma: Optional[dict] = None,
)