# File R/tergm.EGMME.SA.R in package tergm, part of the
# Statnet suite of packages for network analysis, https://statnet.org .
#
# This software is distributed under the GPL-3 license. It is free,
# open source, and has the attribution requirements (GPL Section 7) at
# https://statnet.org/attribution .
#
# Copyright 2008-2023 Statnet Commons
################################################################################
tergm.EGMME.SA <- function(theta0, nw, model, model.mon,
control, proposal, eval.optpars,
verbose=FALSE){
###### Set the constants and convenience variables. ######
offsets <- model$etamap$offsettheta # which parameters are offsets?
p.free <- sum(!model$etamap$offsettheta) # number of free parameters
p <- length(model$etamap$offsettheta) # total number of parameters
p.names <- param_names(model)
q <- length(model.mon$etamap$offsettheta) # number of target statistics
q.names<-param_names(model.mon)
if(control$SA.plot.progress && !dev.interactive(TRUE)){
warning("Progress plot requested on a non-interactive graphics device. Ignoring.")
control$SA.plot.progress <- FALSE # So that we don't print a warning every step.
}
if(verbose) message("Starting optimization with with coef_0 = ( ", paste(theta0, collapse = " "), " ).")
###### Define the optimization run function. ######
do.optimization<-function(states, history, control){
ind.all <- history$ind.all
tid.all <- history$tid.all
oh.all <- history$oh.all
jitters.all <- history$jitters.all
if(verbose) message("Running stochastic optimization... ", appendLF = FALSE)
zs <- if(!is.null(ergm.getCluster(control))){
requireNamespace('parallel')
# Conveniently, the first argument of tergm.EGMME.SA.Phase2.C
# is the state of the optimization, so giving clusterApply a
# list of states will call it for each thread's state.
if(verbose) {message("Calling tergm.EGMME.SA.Phase2.C:"); message_print(gc())}
out <- parallel::clusterApply(ergm.getCluster(control), states, tergm.EGMME.SA.Phase2.C, model, model.mon, proposal, control, verbose=verbose)
if(verbose) message_print(gc())
out
}else{
list(tergm.EGMME.SA.Phase2.C(states[[1]], model, model.mon, proposal, control, verbose=verbose))
}
if(verbose) message("Finished. Extracting.")
for(i in seq_along(states)){
# Extend the observation index and thread id vectors.
ind.all <- c(ind.all, sum(tid.all==i) + 1:(control$SA.runlength*control$SA.interval))
tid.all <- c(tid.all, rep(i, control$SA.runlength*control$SA.interval))
# Extract and store the history of jitters.
jitters.all <- rbind(jitters.all,zs[[i]]$opt.history[,p+1:p,drop=FALSE])
colnames(jitters.all) <- p.names
zs[[i]]$opt.history <- zs[[i]]$opt.history[,-(p+1:p),drop=FALSE]
# Extract and store histroy of trials.
oh.all <- rbind(oh.all,zs[[i]]$opt.history)
colnames(oh.all) <- c(p.names,q.names)
}
# Limit the size of the full history.
# TODO: Store a thinned or summarized version of "forgotten" iterations.
# TODO: Add an option to log it to external storage.
if(max(ind.all)>control$SA.oh.memory){
remember.after.ind <- max(ind.all) - control$SA.oh.memory
# Cut the ind.all last!
tid.all <-tid.all[ind.all>remember.after.ind]
oh.all <- oh.all[ind.all>remember.after.ind,,drop=FALSE]
jitters.all <- jitters.all[ind.all>remember.after.ind,,drop=FALSE]
ind.all <- ind.all[ind.all>remember.after.ind]
# Shift the indices so that they would start at 1.
ind.all <- ind.all - remember.after.ind
}
history$ind.all <- ind.all
history$tid.all <- tid.all
history$jitters.all <- jitters.all
history$oh.all <- oh.all
min.ind.last <- max(ind.all) - control$SA.runlength*control$SA.interval + 1
min.ind.keep <- max(ind.all) - max(
max(ind.all)*control$SA.keep.oh,
min(
max(
control$SA.runlength*control$SA.interval*control$SA.keep.min.runs*length(states),
control$SA.keep.min*length(states)
),
max(ind.all)
)
) + 1
# Extract and store subhistories of interest.
history$ind <- ind <- ind.all[ind.all>=min.ind.keep]
history$tid <- tid <-tid.all[ind.all>=min.ind.keep]
history$ind.last <- ind.last <-ind.all[ind.all>=min.ind.last]
history$tid.last <- tid.last <- tid.all[ind.all>=min.ind.last]
history$oh <- oh <- oh.all[ind.all>=min.ind.keep,,drop=FALSE]
history$oh.last <- oh.last <- oh.all[ind.all>=min.ind.last,,drop=FALSE]
history$jitters <- jitters <- jitters.all[ind.all>=min.ind.keep,,drop=FALSE]
history$jitters.last <- jitters.last <- jitters.all[ind.all>=min.ind.last,,drop=FALSE]
rm(ind.all, tid.all, jitters.all, oh.all); gc()
# Plot if requested.
if(control$SA.plot.progress && dev.interactive(TRUE)){
requireNamespace('lattice')
get.dev("progress.plot")
thin <- (nrow(oh)-1)%/%(control$SA.max.plot.points/length(states)) + 1
cols <- floor(sqrt(ncol(oh)))
layout <- c(cols,ceiling(ncol(oh)/cols))
#' @importFrom coda mcmc.list
suppressWarnings(print(lattice::xyplot(window(do.call(mcmc.list,by(as.data.frame(oh),INDICES=list(tid=tid),mcmc,start=min.ind.keep)), thin=thin), panel = function(...) {lattice::panel.xyplot(...);lattice::panel.abline(0, 0)}, as.table = TRUE, layout = layout, xlab=NULL)))
}
# Extract and return the "states" and the "history".
list(states = lapply(zs, function(z) list(nw = z$newnetwork,
nw.diff = z$nw.diff,
eta = z$eta)
),
history = history)
}
# This function essentially runs the chain forward with gradient etc. set to 0.
do.dummy.run <- function(states, history, control, burnin, steps, eta=NULL){
control$GainM <- matrix(0, nrow=p, ncol=q)
control$dejitter <- matrix(0, nrow=p, ncol=p) # Dejitter tries to cancel the effect of jitter on the optimizer.
control$par.guard <- control$jitter<-rep(0,p)
control$dev.guard <- rep(0,q)
control$dev.guard[] <- 1e10 # A huge value
control$par.guard[!offsets] <- 1e10
control$SA.keep.min.runs <- 0
control$SA.runlength <- 1
control$SA.burnin <- burnin
control$SA.interval <- steps
for(i in seq_along(states)){
if(!is.null(eta)) states[[i]]$eta <- eta
}
do.optimization(states,history,control)
}
interpolate.par <- function(h.fit, w=diag(1,nrow=ncol(h.fit))){
x <- t(h.fit[-1,,drop=FALSE])
y <- -cbind(h.fit[1,])
c(solve(t(x)%*%w%*%x)%*%t(x)%*%w%*%y)
}
#' @importFrom MASS ginv
V.sandwich <- function(w, G, V.stat=ginv(w)){
solve(t(G)%*%w%*%G)%*%t(G)%*%w%*%V.stat%*%w%*%G%*%solve(t(G)%*%w%*%G)
}
best.states <- function(states, history){
w <- ginv(cov(history$oh.all[,-(1:p),drop=FALSE]))
best.i <- which.min(mahalanobis((history$oh.all[-1,-(1:p),drop=FALSE]+history$oh.all[-nrow(history$oh.all),-(1:p),drop=FALSE])/2,0,w,inverted=TRUE))
best.par <- history$oh.all[best.i,1:p][!offsets]
lapply(states, function(state){
if(p.free) state$eta[!model$etamap$offsettheta] <- best.par[seq_len(p.free)]
state
})
}
##### Construct the initial state. ######
history <- list()
history$ind.all <- history$tid.all <- history$oh.all <- history$jitters.all <- state <- NULL
for(restart in 1:control$SA.restarts){
nw.diff <- model.mon$nw.stats - model.mon$target.stats # nw.diff keeps track of the difference between the current network and the target statistics.
states <- replicate(nthreads(control),
{
list(nw=nw,
eta = ergm.eta(theta0, model$etamap),
nw.diff = nw.diff)
},
simplify=FALSE
)
message('======== Phase 1: Burn in, get initial gradient values, and find a configuration under which all targets vary. ========')
###### Set up and run the burn-in. ######
control$changes <- control$collect <- FALSE
control.phase1<-control
control.phase1$time.samplesize <- 1
control.phase1$time.burnin <- control$SA.burnin
control.phase1$time.interval <- 1
message("Burning in... ", appendLF = FALSE)
zs <- if(!is.null(ergm.getCluster(control))){
requireNamespace('parallel')
if(verbose) {message("Calling tergm_MCMC_sample:"); message_print(gc())}
out <- parallel::clusterApply(ergm.getCluster(control), seq_along(states), function(i) tergm_MCMC_sample(states[[i]]$nw, model, model.mon, proposal, eta=states[[i]]$eta, control=control.phase1, verbose=verbose))
if(verbose) message_print(gc())
out
}else{
list(tergm_MCMC_sample(states[[1]]$nw, model, model.mon, proposal, eta=states[[1]]$eta, control=control.phase1, verbose=verbose))
}
message("Done.")
# Update the state with burn-in results.
for(i in seq_along(zs)){
states[[i]]$nw <- zs[[i]]$newnetwork
states[[i]]$nw.diff <- states[[i]]$nw.diff + zs[[i]]$statsmatrix.mon[NROW(zs[[i]]$statsmatrix.mon),]
}
###### Gradient estimation and getting to a decent configuration. ######
## Here, "decent configuration" is defined as one where all target
## statistics have some variability --- none are stuck.
# Start with jitter-only.
control$GainM <- matrix(0, nrow=p, ncol=q)
control$dejitter <- matrix(0, nrow=p, ncol=p) # Dejitter tries to cancel the effect of jitter on the optimizer.
control$par.guard <- control$jitter<-rep(0,p)
control$dev.guard <- rep(0,q)
control$jitter[!offsets] <- control$SA.phase1.jitter
control$dev.guard[] <- 1e10 # A huge value
control$par.guard[!offsets] <- control$SA.phase1.jitter * 4
## Adjust the number of time steps between jumps using burn-in.
edge.ages <- unlist(sapply(states, function(state) (((if(is(state$nw, "ergm_state")) state$nw$nw0 else state$nw) %n% "time") - edgelist_with_lasttoggle(state$nw)[,3] + 1)))
control$SA.burnin<-control$SA.interval<- round(min(control$SA.max.interval, max(control$SA.min.interval, if(length(edge.ages)>0) control$SA.interval.mul*mean(edge.ages)))/2)
if(is.nan(control$SA.burnin)|is.null(control$SA.burnin)|is.na(control$SA.burnin))
control$SA.burnin <- control$SA.interval <- 10 # TODO: Kirk : check this
if(verbose>1){
message("New interval: ", control$SA.interval)
}
for(try in 1:control$SA.phase1.tries){
if(verbose) message('======== Attempt ', try, ' ========') else message('Attempt ', try,' :')
for(run in 1:control$SA.phase1.minruns){
tmp <- if(control$SA.restart.on.err) try(do.optimization(states, history, control), silent=!verbose) else do.optimization(states, history, control)
if(inherits(tmp, "try-error") || all(apply(tmp$history$oh.last[,-(1:p),drop=FALSE],2,var)<sqrt(.Machine$double.eps))){
message("Something went very wrong. Restarting with smaller gain.")
control$SA.init.gain <- control$SA.init.gain * control$SA.gain.decay
do.restart <- TRUE
break
}else{
states <- tmp$states
history <- tmp$history
do.restart <- FALSE
}
}
states <- best.states(states, history)
control$gain <- control$SA.init.gain
out <- if(control$SA.restart.on.err) try(eval.optpars(states, history, control, TRUE,restart>1,FALSE), silent=!verbose) else eval.optpars(states, history, control, TRUE,restart>1,FALSE)
if(inherits(out, "try-error") || all(apply(history$oh.last[,-(1:p),drop=FALSE],2,var)<sqrt(.Machine$double.eps))){
message("Something went very wrong. Restarting with smaller gain.")
control$SA.init.gain <- control$SA.init.gain * control$SA.gain.decay
do.restart <- TRUE
break
}else do.restart <- FALSE
control <- out$control
if(all(!out$ineffectual.pars) && all(!out$bad.fits)){
message("All parameters have some effect and all statistics are moving. Proceeding to Phase 2.")
break
}
if(try==control$SA.phase1.tries) stop("The optimizer was unable to find a reasonable configuration: one or more statistics are still stuck after multiple tries, and one or more parameters do not appear to have any robust effect.")
}
if(do.restart) next
###### Main optimization run. ######
message('======== Phase 2: Find and refine the estimate. ========')
for(subphase in 1:control$SA.phase2.levels.max){
if(verbose) message('======== Subphase ', subphase, ' ========') else message('Subphase 2.', subphase, ' ', appendLF=FALSE)
control$gain <- control$SA.init.gain*control$SA.gain.decay^(subphase-1)
stepdown.count <- control$SA.stepdown.ct
for(regain in 1:control$SA.phase2.repeats){
tmp <- if(control$SA.restart.on.err) try(do.optimization(states, history, control), silent=!verbose) else do.optimization(states, history, control)
if(inherits(tmp, "try-error") || all(apply(tmp$history$oh.last[,-(1:p),drop=FALSE],2,var)<sqrt(.Machine$double.eps))){
message("Something went very wrong. Restarting with smaller gain.")
control$SA.init.gain <- control$SA.init.gain * control$SA.gain.decay
do.restart <- TRUE
break
}else{
states <- tmp$states
history <- tmp$history
do.restart <- FALSE
}
if(verbose){
message("New parameters:")
for(state in states) message_print(state$eta)
}
## Get updated gain and other values
out <- if(control$SA.restart.on.err) try(eval.optpars(states, history, control, FALSE,TRUE,TRUE), silent=!verbose) else eval.optpars(states, history, control, FALSE,TRUE,TRUE)
if(inherits(out, "try-error") || all(apply(history$oh.last[,-(1:p),drop=FALSE],2,var)<sqrt(.Machine$double.eps))){
message("Something went very wrong. Restarting with smaller gain.")
control$SA.init.gain <- control$SA.init.gain * control$SA.gain.decay
do.restart <- TRUE
break
}else do.restart <- FALSE
control <- out$control
## Run two tests here:
## 1) Are the estimating equations, on average, 0?
## 2) Does there appear to be a trend in their sum of squares?
# Calculate the approximate estimating equation values:
ys <- history$oh[,-(1:p),drop=FALSE]%*%out$w%*%out$G
# This test is fast, so no need to thin.
p.val.1 <- try(approx.hotelling.diff.test(ys)$p.value)
if(is.na(p.val.1)) p.val.1 <- 0
# Thin the data to keep from bogging down.
x <- unique(round(seq(from=1,to=NROW(history$oh),length.out=control$SA.stepdown.maxn)))
y <- sqrt(mahalanobis(ys,0,ginv(cov(ys)),inverted=TRUE)[x])
i <- history$ind[x]
t <- history$tid[x]
for(thread in unique(t)) i[t==thread] <- rank(i[t==thread])
#' @importFrom nlme gls corAR1
fit.2 <- try(summary(gls(y~x,correlation=corAR1(form=~i|t)))$tTable[2,c(1,4)])
if(!inherits(p.val.1, "try-error") && !inherits(fit.2, "try-error")){
p.val.2 <- fit.2[2]/2 # We are interested in one-sided decline here.
est.2 <- fit.2[1]
if(est.2>0) p.val.2 <- 1-p.val.2 # If it's actually getting worse one-sided p-value is thus.
if(verbose){
message("Estimating equations = 0 p-value: ", p.val.1, " , trending: ", p.val.2, " .")
}
p.vals <- c(p.val.1,p.val.2)
fisher.pval <- function(p.vals){
p.vals <- unlist(p.vals)
df <- 2*length(p.vals)
pchisq(-2*sum(log(p.vals)), df, lower.tail=FALSE)
}
if(fisher.pval(p.vals)>control$SA.stepdown.p){
stepdown.count <- stepdown.count - 1
if(stepdown.count<=0){
if(verbose) message("Estimating equations do not significantly differ from 0 and neither they nor the parameters exhibit a significant trend. Reducing gain.")
else message("\\", appendLF = FALSE)
stepdown.count <- control$SA.stepdown.ct
if(!verbose) message("")
break
}else{
if(verbose) message("Estimating equations do not significantly differ from 0 and do not exhibit a significant trend. ", stepdown.count, " / ", control$SA.stepdown.ct, " to go.")
else message("\\", appendLF = FALSE)
}
}else{
stepdown.count <- control$SA.stepdown.ct
if(verbose) message("Estimating equations significantly differ from 0 or exhibit a significant trend. Resetting counter.")
else message("/", appendLF = FALSE)
}
}else{
if(verbose) message("Problem testing estimating equations. Continuing with current gain.")
else message("!/", appendLF = FALSE)
}
if(do.restart) break
}
if(do.restart) break
#### Decide whether to stop.
## Run through minimal number of phases.
if(subphase<control$SA.phase2.levels.min) next
## Run three tests here:
## 1) Is the stochastic approximation estimate of the GMM sufficiently precise?
## 2) Is there strong evidence of nonlinearity?
## 3) Does there appear to be a trend in parameter values?
# Test 1:
par.se <- {
h <- history$oh[,1:p,drop=FALSE][,!offsets,drop=FALSE]
apply(h,2,function(x) sd(x)/sqrt(effectiveSize(x)))
}
sandwich.se <- sqrt(diag(V.sandwich(out$w,out$G,out$v)))
if(verbose){
message("Approximate standard error of the estimate:")
message_print(sandwich.se)
message("Approximate standard error of window means:")
message_print(par.se)
message("par. var. / (std. var. + par. var.):")
message_print(par.se^2/(sandwich.se^2+par.se^2))
}
# Test 2:
ys <- history$oh[, -(1:p), drop=FALSE]
xs <- history$oh[, 1:p, drop=FALSE][,!offsets,drop=FALSE]
nlin.totest <- rep(c(FALSE,TRUE),c(p.free+1, p.free*2))
nlin.fit <- lm(ys~xs+I(xs^2)+I(xs^3))
nlin.nfs <- matrix(apply(cbind(resid(nlin.fit)),2,function(x) sum(tapply(x,list(history$tid),length))/sum(tapply(x,list(history$tid),effectiveSize))), nrow=1+p.free*3, ncol=q, byrow=TRUE)
nlin.coef <- cbind(coef(nlin.fit))
nlin.vcov <- vcov(nlin.fit)
drop <- apply(is.na(nlin.coef),1,any)
nlin.coef <- nlin.coef[nlin.totest & !drop,,drop=FALSE]
nlin.nfs <- nlin.nfs[nlin.totest & !drop,,drop=FALSE]
nlin.vcov <- t(nlin.vcov[nlin.totest & !drop, nlin.totest & !drop, drop=FALSE]*sqrt(c(nlin.nfs)))*sqrt(c(nlin.nfs))
chi2 <- mahalanobis(c(nlin.coef),0,ginv(nlin.vcov),inverted=TRUE)
p.val.1 <- pchisq(chi2, length(nlin.coef), lower.tail=FALSE)
if(verbose) message("Local nonlinearity p-value: ", p.val.1)
if(subphase == control$SA.phase2.levels.max){
if(verbose) message("Maximum number of gain levels exceeded. Stopping.", appendLF = FALSE)
}else{
if(all(par.se^2/(sandwich.se^2+par.se^2) > control$SA.phase2.max.mc.se)){
if(verbose) message("EGMME does not appear to be estimated with sufficient prescision. Continuing.")
next
}
if(p.val.1 < control$SA.stop.p){
if(verbose) message("There is evidence of local nonlinearity. Continuing.")
next
}
}
###### If we've gotten this far, proceed to Part 3. ######
## Refine the estimate.
if(inherits(state, "try-error")) stop("Something went wrong too many times. Try better starting values or reducing control$SA.init.gain.")
eta <- states[[1]]$eta
eta.free <- switch(control$SA.refine,
mean = colMeans(history$oh[,1:p,drop=FALSE][,!offsets,drop=FALSE]),
linear = interpolate.par(out$oh.fit,out$w),
none = if(is.null(ergm.getCluster(control))) eta[!offsets] else stop("No interpolation does not make sense with multithreaded fitting."))
if(p.free) eta[!model$etamap$offsettheta] <- eta.free[seq_len(p.free)]
if(verbose){
message("Refining the estimate using the ", control$SA.refine, " method. New estimate:")
message_print(eta)
}
## Sample to estimate standard error and test if we've actually arrived.
G <- out$G
w <- out$w
V.par<-matrix(NA,p,p)
if(control$SA.se){
message('======== Phase 3: Simulate from the fit and estimate standard errors. ========')
if(verbose) message("Evaluating target statistics at the estimate.")
# This is to avoid the latest sample "pushing out" everything else.
control$SA.keep.min <- round(nrow(history$oh)/length(states))+control$SA.phase3.samplesize.runs*control$SA.runlength*control$SA.interval
tmp <- do.dummy.run(states, history, control, control$SA.burnin, control$SA.phase3.samplesize.runs*control$SA.runlength*control$SA.interval, eta=eta)
if(verbose) message("Finished.")
states <- tmp$states
history <- tmp$history
rm(tmp); gc()
min.ind.se <- max(history$ind.all)-control$SA.phase3.samplesize.runs*control$SA.runlength*control$SA.interval + 1
sm.mon <- history$oh.all[history$ind.all>=min.ind.se,-(1:p),drop=FALSE]
sm.tid <- history$tid.all[history$ind.all>=min.ind.se]
V.stat<-cov(sm.mon)
ee <- sm.mon %*% ginv(V.stat) %*% G
p.val.1 <- approx.hotelling.diff.test(ee)$p.value
if(verbose) message("Estimating equation = 0 p-value: ", p.val.1)
if(subphase == control$SA.phase2.levels.max){
if(verbose) message("Maximum number of gain levels exceeded. Stopping.", appendLF=FALSE)
}else{
if(p.val.1 < control$SA.stop.p){
if(verbose) message("Simulated values of estimating equations are not centered around 0. Continuing.")
next
}else{
if(verbose) message("Simulated values of estimating equations are centered around 0. Stopping.")
}
}
V.par[!offsets,!offsets]<-V.sandwich(w,G,V.stat)
sm.mons <- lapply(unique(sm.tid), function(t) sm.mon[sm.tid==t,,drop=FALSE])
sm.mon <- do.call(mcmc.list, lapply(sm.mons, mcmc, start=control$SA.burnin))
}else{
V.par[!offsets,!offsets]<-V.sandwich(w,G,out$v)
sm.mon <- NULL
}
mc.se <- rep(NA, p)
mc.se[!offsets] <- par.se
break
}
if(!do.restart) break # If we've gotten this far, no restart conditions have been triggered, so we are good.
}
list(newnetwork=if(control$SA.se) zs[[1]]$newnetwork else states[[1]]$nw,
newnetworks=if(control$SA.se) lapply(zs,"[[","newnetwork") else lapply(states,"[[","nw"),
init=theta0,
covar=V.par,
mc.se = mc.se,
eta=eta,
opt.history=history$oh.all,
sample=sm.mon,
network=nw)
}
tergm.EGMME.SA.Phase2.C <- function(state, model, model.mon,
proposal, control, verbose) {
on.exit(ergm_Cstate_clear())
model.comb <- c(model, model.mon)
proposal$aux.slots <- model.comb$slots.extra.aux$proposal
eta.comb <- c(deInf(state$eta), rep(0,model.mon$etamap$etalength))
ergmstate <- ergm_state(state$nw, model=model.comb, proposal=proposal)
maxedges <- max(NVL(control$MCMC.maxedges, Inf), network.edgecount(ergmstate))
maxchanges <- max(control$MCMC.maxchanges, network.edgecount(ergmstate))
z <- .Call("MCMCDynSArun_wrapper",
ergmstate,
as.integer(nparam(model.mon,canonical=TRUE)),
# Parameter fitting.
as.double(eta.comb),
as.double(state$nw.diff),
as.integer(control$SA.runlength),
as.double(control$GainM),
as.double(control$jitter),
as.double(control$dejitter), # Add a little bit of noise to parameter guesses.
as.double(control$dev.guard),
as.double(control$par.guard),
# MCMC settings.
as.integer(control$SA.burnin),
as.integer(control$SA.interval),
as.integer(control$MCMC.burnin.min),
as.integer(control$MCMC.burnin.max),
as.double(control$MCMC.burnin.pval),
as.double(control$MCMC.burnin.add),
as.integer(deInf(maxedges, "maxint")),
as.integer(maxchanges),
as.integer(max(verbose-1,0)),
PACKAGE="tergm")
if(z$status != 0) stop("DynSA errored with code ", z$status)
z$state <- update(z$state)
eta <- z$eta[seq_len(model$etamap$etalength)]
names(eta) <- param_names(model, canonical = TRUE)
newnetwork <- as.network(z$state)
list(nw.diff=z$nw.diff,
newnetwork=newnetwork,
eta=eta,
opt.history=matrix(z$opt.history,ncol=2*model.comb$etamap$etalength - model.mon$etamap$etalength,byrow=TRUE))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.