---
title: "Introduction to cusna"
output: litedown::html_format
vignette: >
  %\VignetteEngine{litedown::vignette}
  %\VignetteIndexEntry{Introduction to cusna}
  %\VignetteEncoding{UTF-8}
---

`cusna` gives R a native, bit-reproducible engine for stochastic
actor-oriented models (SAOM --- the model family of
[RSiena](https://cran.r-project.org/package=RSiena)) and exponential random
graph models. The engine is a C ABI over CUDA kernels and C++ host logic,
compiled with the package --- **no Python runtime**. Following the design of
the [`torch`](https://cran.r-project.org/package=torch) package, the CRAN
build is CPU-only from source; the GPU path is compiled from source when a
CUDA 12.x toolkit is detected at configure time.

Because the engine ships inside the package, the estimation chunks in this
vignette actually run (on small synthetic data, with reduced simulation
counts to keep the build fast).

## Setup

Nothing to set up --- the engine is already compiled in:

```{r}
library(cusna)
cusna_abi_version()
cusna_has_cuda()      # FALSE on the CPU-only build

# Keep the vignette build well within CRAN's two-core policy: the native CPU
# backend is OpenMP-parallel over simulation chains, so run it single-threaded
# here. cusna_set_threads() calls the OpenMP runtime directly, so it takes
# effect immediately (unlike Sys.setenv(OMP_NUM_THREADS = ...), whose effect
# on an already-initialised OpenMP runtime is platform-dependent). The demo
# data is tiny, so a single thread is still fast.
old_threads <- cusna_set_threads(1L)
```

## Specifying a model

Model terms are plain R objects:

```{r}
effects <- list(
  cusna_effect("density"),
  cusna_effect("recip"),
  cusna_effect("transTrip")
)
effects[[1]]
```

Covariate effects take a covariate; the engine centers it for you. A covariate
is a numeric vector of length *n* (constant) or an *n*-by-waves matrix
(changing):

```{r}
alcohol <- c(3, 2, 2, 1, 3)         # one value per actor (toy)
cusna_effect("egoX", covariate = alcohol)
```

User interactions multiply the change statistics of two or three component
effects, which must appear earlier in the list. Component positions are
1-based:

```{r}
list(
  cusna_effect("egoX", covariate = alcohol),
  cusna_effect("altX", covariate = alcohol),
  cusna_interaction(c(1, 2))          # egoX x altX
)
```

## Data and estimation

Build a panel from a list of adjacency matrices (missing ties are `NA` or the
code 9; structural zeros/ones are 10/11). `saom_data()` applies RSiena's
imputation and per-period masks --- a construction validated bit-for-bit
against the reference implementation:

```{r}
set.seed(7)
n <- 20
w1 <- matrix(as.integer(runif(n * n) < 0.12), n, n); diag(w1) <- 0L
w2 <- w1; flip <- sample(n * n, 40); w2[flip] <- 1L - w2[flip]; diag(w2) <- 0L
w3 <- w2; flip <- sample(n * n, 40); w3[flip] <- 1L - w3[flip]; diag(w3) <- 0L

dat <- saom_data(list(w1, w2, w3))
dat
```

Estimation runs the RSiena-style Robbins--Monro method of moments on the
native simulator. (The tiny `mom_control()` counts below keep this vignette
quick; the defaults follow RSiena and are the ones to use in practice.)

```{r}
fit <- mom_estimate(
  dat,
  effects = list(cusna_effect("density"), cusna_effect("recip")),
  conditional = TRUE,                 # RSiena's default for one network
  control = mom_control(n1 = 100, nsub = 1, n2 = 10, batch2 = 50, n3 = 200))
fit
```

The result is an ordinary R object:

```{r}
summary(fit)          # RSiena-style table (rates, estimates, s.e., t-conv)
coef(fit)             # named parameter vector
round(vcov(fit), 4)   # delta-method covariance
```

### Covariates and behavior co-evolution

Actor covariates enter through the covariate effects; network--behavior
co-evolution adds a `behavior` matrix and behavior effects, and uses
unconditional estimation (as in RSiena for more than one dependent variable):

```{r, eval=FALSE}
dat <- saom_data(list(w1, w2, w3), behavior = drink)   # drink: n-by-waves

fit <- mom_estimate(
  dat,
  effects = list(cusna_effect("density"), cusna_effect("recip"),
                 cusna_effect("transTrip"),
                 cusna_effect("altX", dyn = TRUE),      # covariate = the behavior
                 cusna_effect("egoX", dyn = TRUE)),
  beh_effects = list(cusna_beh_effect("linear"), cusna_beh_effect("quad"),
                     cusna_beh_effect("avSim")),
  conditional = FALSE)
```

Several dependent networks over one actor set co-evolve through
`saom_multinet_data()` + `mom_estimate_multinet()`, with cross-network
effects (`"crprod"`, `"from"`, ...) pointing at their source network via
`net_ref`.

## Tuning the estimator

`mom_control()` exposes the Robbins--Monro schedule. The defaults follow
RSiena; raise `diagonalize` towards 1 and add phase-3 simulations for harder,
nearly collinear specifications:

```{r}
str(mom_control())
```

| Argument | Meaning |
|---|---|
| `firstg` | initial gain |
| `nsub`, `n2` | phase-2 subphases and their iteration counts |
| `batch2` | simulation batch per phase-2 iteration |
| `n1`, `nD`, `n3` | phase-1 / derivative / phase-3 simulation counts |
| `diagonalize` | shrink the derivative towards its diagonal (0--1) |
| `seed` | random seed |

On nearly collinear pairs (e.g. `transTrip` + `gwesp`) plain Robbins--Monro
may not converge from a cold start: increase `diagonalize`, or drive RSiena's
own estimator through the native simulator with `cusna_fran()` --- see
`vignette("siena07-backend")`.

## Beyond SAOM: ERGM, TERGM, STERGM, ALAAM

The same engine serves cross-sectional and temporal exponential random graph
models, and autologistic actor attribute models. Sufficient statistics and a
TNT sampler are direct calls:

```{r}
x <- w1                                       # any 0/1 adjacency
terms <- list(ergm_term("edges"), ergm_term("mutual"))
ergm_stats(x, terms, directed = TRUE)

ergm_simulate(x, coef = c(-2, 1), terms = terms, nsim = 5,
              directed = TRUE, seed = 1)
```

Estimators (each validated against its R reference):

```{r, eval=FALSE}
ergm_mcmle(x, terms, directed = TRUE)       # MCMC-MLE       (matches ergm::ergm)
ergm_mple(x)                                # pseudo-likelihood demo
tergm_mple(list(w1, w2, w3))                # TERGM           (matches btergm)
tergm_simulate(w2, lag = w1, coef = c(-4, 3, 2, 0.5))
stergm_cmle(list(w1, w2, w3),               # separable TERGM (matches tergm CMLE)
            formation = list(ergm_term("edges"), ergm_term("mutual")))
alaam_mple(y, net = x)                      # ALAAM MPLE      (matches glm)
alaam_mcmle(y, net = x)                     # ALAAM MCMC-MLE
```

## Validation

The native stack is validated in layers: the CUDA kernels are bit-identical
to the reference GPU implementation; the host statistics reproduce RSiena
moment targets on public datasets (s50, Knecht, the COW alliance panel) to
machine precision; the R-side data preparation, effect preprocessing, and
moment targets are bit-identical to the reference on seven datasets; and the
full Robbins--Monro estimator agrees with the reference estimator within
simulation standard errors on baseline, covariate, co-evolution, and
multi-network models. The ERGM samplers match `ergm::simulate`
distributionally, and the ERGM/TERGM/STERGM/ALAAM estimators match
`ergm`/`btergm`/`tergm`/`glm` on their benchmark models.

An NVIDIA GPU is optional throughout: the CPU build compiles from source on
any C++17 toolchain, and `cusna_has_cuda()` reports which build you have.