Skip to content

Commit f316d81

Browse files
krivitkrivit
committed
Undid whitespace changes in the pull request.
1 parent e388ec6 commit f316d81

File tree

2 files changed

+56
-57
lines changed

2 files changed

+56
-57
lines changed

‎R/tergm.getMCMCsample.R

Lines changed: 18 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -9,15 +9,15 @@
99
################################################################################
1010

1111
#' Collects a sample of networks and returns the statistics of each sample
12-
#'
12+
#'
1313
#' \code{tergm_MCMC_sample} is a low-level internal function not intended to
1414
#' be called directly by end users. It collects a sample of networks and
15-
#' returns the statistics of each sample, along with a toggle matrix of the
15+
#' returns the statistics of each sample, along with a toggle matrix of the
1616
#' changes needed from the original network to each in the sample.
17-
#'
17+
#'
1818
#' This function is normally called inside [simulate.tergm()] functions
1919
#' to prepare inputs for the C sampling code and return its results
20-
#'
20+
#'
2121
#' @aliases tergm_MCMC_sample tergm_MCMC_slave
2222
#' @param nw a [`network`] object
2323
#' @param model the model, as returned by [`ergm_model`]
@@ -27,14 +27,14 @@
2727
#' @param eta the vector of natural parameters
2828
#' @param control the list of control parameters
2929
#' @template verbose
30-
#' @return returns the MCMC sample as a list containing:
30+
#' @return returns the MCMC sample as a list containing:
3131
#' \itemize{
3232
#' \item statsmatrix.gen: the matrix of sampled statistics for \code{model},
3333
#' relative to the initial network
3434
#' \item statsmatrix.mon: the matrix of sampled statistics for \code{model.mon},
35-
#' relative to the initial network
35+
#' relative to the initial network
3636
#' \item newnetwork: \code{ergm_state} with the final network from the
37-
#' sampling process
37+
#' sampling process
3838
#' \item changed: a matrix of changes, where the first column is
3939
#' the timestamp of the change, the second and third columns are the tail and head
4040
#' (respectively) of the changed dyad, and the fourth column is the edge state to which
@@ -50,22 +50,22 @@ tergm_MCMC_sample <- function(nw, model, model.mon = NULL,
5050
verbose=FALSE,...,
5151
eta = ergm.eta(theta, model$etamap)
5252
){
53-
# this is where we combine models and pad out eta
53+
# this is where we combine models and pad out eta
5454
# with 0s as necessary to accomodate the monitoring model
5555
model.comb <- c(model, model.mon)
5656
proposal$aux.slots <- model.comb$slots.extra.aux$proposal
57-
57+
5858
eta.comb <- c(eta, rep(0, NVL(model.mon$etamap$etalength, 0)))
5959

6060
# always collect if monitoring model is passed
6161
control$collect <- NVL(control$collect, TRUE) || !is.null(model.mon)
6262

6363
#
6464
# Check for truncation of the returned edge list
65-
#
66-
65+
#
66+
6767
state <- ergm_state(nw, model=model.comb, proposal=proposal, stats=rep(0,nparam(model.comb, canonical=TRUE)))
68-
68+
6969
z <- tergm_MCMC_slave(state, eta.comb, control, verbose)
7070

7171
if(z$status)
@@ -77,7 +77,7 @@ tergm_MCMC_sample <- function(nw, model, model.mon = NULL,
7777
)
7878

7979
state <- z$state
80-
80+
8181
diffedgelist<-if(control$changes) {
8282
if(z$diffnwtime[1]>0){
8383
tmp <- cbind(z$diffnwtime[2:(z$diffnwtime[1]+1)],z$diffnwtails[2:(z$diffnwtails[1]+1)],z$diffnwheads[2:(z$diffnwheads[1]+1)],z$diffnwdirs[2:(z$diffnwdirs[1]+1)])
@@ -94,7 +94,7 @@ tergm_MCMC_sample <- function(nw, model, model.mon = NULL,
9494
mode(diffedgelist) <- "integer" # Might save some memory.
9595

9696
statsmatrix <- z$statsmatrix
97-
97+
9898
if(!is.null(statsmatrix)) colnames(statsmatrix) <- param_names(model.comb, canonical = TRUE)
9999

100100
# this is where we separate monitored stats from generative stats if model.mon is passed
@@ -125,7 +125,7 @@ tergm_MCMC_slave <- function(state, eta, control, verbose){
125125

126126
maxedges <- NVL(control$MCMC.maxedges, Inf)
127127
maxchanges <- control$MCMC.maxchanges
128-
128+
129129
z <- .Call("MCMCDyn_wrapper",
130130
state,
131131
as.double(deInf(eta)),
@@ -146,16 +146,16 @@ tergm_MCMC_slave <- function(state, eta, control, verbose){
146146
PACKAGE="tergm")
147147

148148
if(z$status) return(z) # If there is an error.
149-
149+
150150
z$state <- update(z$state)
151-
151+
152152
statsmatrix <-
153153
if(collect) matrix(z$s, nrow=control$time.samplesize+1,
154154
ncol=nparam(state,canonical=TRUE),
155155
byrow = TRUE)[-1,,drop=FALSE]
156156
else
157157
NULL
158-
158+
159159
c(z,
160160
list(statsmatrix = statsmatrix))
161161
}

‎src/MCMCDyn.c

Lines changed: 38 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ SEXP MCMCDyn_wrapper(SEXP stateR, // ergm_state
5151
SEXP sample = PROTECT(allocVector(REALSXP, (asInteger(nsteps) + 1)*m->n_stats));
5252
memset(REAL(sample), 0, (asInteger(nsteps) + 1)*m->n_stats*sizeof(double));
5353
memcpy(REAL(sample), s->stats, m->n_stats*sizeof(double));
54-
54+
5555
kvint difftime, difftail, diffhead, diffto;
5656
kv_init(difftime);
5757
kv_init(difftail);
@@ -64,7 +64,6 @@ SEXP MCMCDyn_wrapper(SEXP stateR, // ergm_state
6464
kv_push(int, diffhead, 0);
6565
kv_push(int, diffto, 0);
6666

67-
6867
SEXP status;
6968
if(MHp) status = PROTECT(ScalarInteger(MCMCSampleDyn(s,
7069
dur_inf,
@@ -73,18 +72,18 @@ SEXP MCMCDyn_wrapper(SEXP stateR, // ergm_state
7372
asInteger(nsteps), asInteger(min_MH_interval), asInteger(max_MH_interval), asReal(MH_pval), asReal(MH_interval_add), asInteger(burnin), asInteger(interval),
7473
asInteger(verbose))));
7574
else status = PROTECT(ScalarInteger(MCMCDyn_MH_FAILED));
76-
75+
7776
const char *outn[] = {"status", "s", "state", "diffnwtime", "diffnwtails", "diffnwheads", "diffnwdirs", ""};
7877
SEXP outl = PROTECT(mkNamed(VECSXP, outn));
7978
SET_VECTOR_ELT(outl, 0, status);
8079
SET_VECTOR_ELT(outl, 1, sample);
81-
80+
8281
/* record new generated network to pass back to R */
8382
if(asInteger(status) == MCMCDyn_OK){
8483
s->stats = REAL(sample) + asInteger(nsteps)*m->n_stats;
8584
SET_VECTOR_ELT(outl, 2, ErgmStateRSave(s));
8685
}
87-
86+
8887
SET_VECTOR_ELT(outl, 3, PROTECT(kvint_to_SEXP(difftime)));
8988
SET_VECTOR_ELT(outl, 4, PROTECT(kvint_to_SEXP(difftail)));
9089
SET_VECTOR_ELT(outl, 5, PROTECT(kvint_to_SEXP(diffhead)));
@@ -95,7 +94,7 @@ SEXP MCMCDyn_wrapper(SEXP stateR, // ergm_state
9594
kv_destroy(diffhead);
9695
kv_destroy(diffto);
9796

98-
ErgmStateDestroy(s);
97+
ErgmStateDestroy(s);
9998
PutRNGstate(); /* Disable RNG before returning */
10099
UNPROTECT(7);
101100
return outl;
@@ -107,9 +106,9 @@ SEXP MCMCDyn_wrapper(SEXP stateR, // ergm_state
107106
Using the parameters contained in the array eta, obtain the
108107
network statistics for a sample of size nsteps. burnin is the
109108
initial number of Markov chain steps before sampling anything
110-
and interval is the number of MC steps between successive
109+
and interval is the number of MC steps between successive
111110
networks in the sample. Put all the sampled statistics into
112-
the statistics array.
111+
the statistics array.
113112
*********************/
114113
MCMCDynStatus MCMCSampleDyn(ErgmState *s,
115114
StoreTimeAndLasttoggle *dur_inf,
@@ -122,7 +121,7 @@ MCMCDynStatus MCMCSampleDyn(ErgmState *s,
122121
kvint *difftime, kvint *difftail, kvint *diffhead, kvint *diffto,
123122
// MCMC settings.
124123
unsigned int nsteps, unsigned int min_MH_interval, unsigned int max_MH_interval, double MH_pval, double MH_interval_add,
125-
unsigned int burnin, unsigned int interval,
124+
unsigned int burnin, unsigned int interval,
126125
// Verbosity.
127126
int verbose){
128127
Network *nwp = s->nwp;
@@ -135,8 +134,8 @@ MCMCDynStatus MCMCSampleDyn(ErgmState *s,
135134
/*if (verbose)
136135
Rprintf("Total m->n_stats is %i; total nsteps is %d\n",
137136
m->n_stats,nsteps);*/
138-
139-
137+
138+
140139
/* Burn in step. */
141140

142141
for(i=0;i<burnin;i++){
@@ -150,18 +149,18 @@ MCMCDynStatus MCMCSampleDyn(ErgmState *s,
150149
// Check that we didn't run out of log space.
151150
if(status==MCMCDyn_TOO_MANY_CHANGES)
152151
return MCMCDyn_TOO_MANY_CHANGES;
153-
152+
154153
// If we need to return a network, then stop right there, since the network is too big to return, so stop early.
155154
if(maxedges!=0 && EDGECOUNT(nwp) >= maxedges-1)
156155
return MCMCDyn_TOO_MANY_EDGES;
157156
}
158-
157+
159158
//Rprintf("MCMCSampleDyn post burnin numdissolve %d\n", *numdissolve);
160-
159+
161160
if (verbose){
162161
Rprintf("Returned from STERGM burnin\n");
163162
}
164-
163+
165164
/* Now sample networks */
166165
for (i=0; i < nsteps; i++){
167166
/* Set current vector of stats equal to previous vector */
@@ -181,16 +180,16 @@ MCMCDynStatus MCMCSampleDyn(ErgmState *s,
181180
stats,
182181
maxchanges, log_changes ? &nextdiffedge : NULL, difftime, difftail, diffhead, diffto,
183182
min_MH_interval, max_MH_interval, MH_pval, MH_interval_add, verbose);
184-
183+
185184
// Check that we didn't run out of log space.
186185
if(status==MCMCDyn_TOO_MANY_CHANGES)
187186
return MCMCDyn_TOO_MANY_CHANGES;
188-
187+
189188
// If we need to return a network, then stop right there, since the network is too big to return, so stop early.
190189
if(maxedges!=0 && EDGECOUNT(nwp) >= maxedges-1)
191190
return MCMCDyn_TOO_MANY_EDGES;
192191
}
193-
192+
194193
//Rprintf("MCMCSampleDyn loop numdissolve %d\n", *numdissolve);
195194
if (verbose){
196195
if( ((3*i) % nsteps)<3 && nsteps > 500){
@@ -241,7 +240,7 @@ MCMCDynStatus MCMCDyn1Step(ErgmState *s,
241240
// MCMC settings.
242241
unsigned int min_MH_interval, unsigned int max_MH_interval, double MH_pval, double MH_interval_add,
243242
// Verbosity.
244-
int verbose){
243+
int verbose){
245244
Network *nwp = s->nwp;
246245
Model *m = s->m;
247246
MHProposal *MHp = s->MHp;
@@ -254,20 +253,20 @@ MCMCDynStatus MCMCDyn1Step(ErgmState *s,
254253
if(stats) addonto(stats, m->workspace, m->n_stats);
255254

256255
/* Run the process. */
257-
256+
258257
double cutoff;
259-
double
258+
double
260259
si = 0, // sum of increments
261260
si2 = 0, // sum of squared increments
262-
sw = 0, // sum of weights
261+
sw = 0, // sum of weights
263262
sw2 = 0 // sum of squared weights
264263
;
265264
double sdecay = 1 - 1.0/min_MH_interval;
266-
265+
267266
unsigned int step=0; // So that we could print out the number of steps later.
268267
for(unsigned int finished = 0, extrasteps = 0; step < max_MH_interval && finished < extrasteps+1; step++) {
269268
unsigned int prev_discord = kh_size(discord);
270-
269+
271270
MHp->logratio = 0;
272271
(*(MHp->p_func))(MHp, nwp); /* Call MHp function to propose toggles */
273272
// Rprintf("Back from proposal; step=%d\n",step);
@@ -277,32 +276,32 @@ MCMCDynStatus MCMCDyn1Step(ErgmState *s,
277276
switch(MHp->togglehead[0]){
278277
case MH_UNRECOVERABLE:
279278
error("Something very bad happened during proposal. Memory has not been deallocated, so restart R soon.");
280-
279+
281280
case MH_IMPOSSIBLE:
282281
Rprintf("MH MHProposal function encountered a configuration from which no toggle(s) can be proposed.\n");
283282
return MCMCDyn_MH_FAILED;
284-
283+
285284
case MH_UNSUCCESSFUL:
286285
case MH_CONSTRAINT:
287286
continue;
288287
}
289288
}
290289

291290
ChangeStats(MHp->ntoggles, MHp->toggletail, MHp->togglehead, nwp, m);
292-
293-
// Rprintf("change stats:");
291+
292+
// Rprintf("change stats:");
294293
/* Calculate inner product */
295294
double ip = dotprod(eta, m->workspace, m->n_stats);
296-
// Rprintf("%f ", m->workspace[i]);
295+
// Rprintf("%f ", m->workspace[i]);
297296
//}
298-
// Rprintf("\n ip %f dedges %f\n", ip, m->workspace[0]);
297+
// Rprintf("\n ip %f dedges %f\n", ip, m->workspace[0]);
299298
/* The logic is to set exp(cutoff) = exp(ip) * qratio ,
300299
then let the MHp probability equal min{exp(cutoff), 1.0}.
301300
But we'll do it in log space instead. */
302301
cutoff = ip + MHp->logratio;
303-
302+
304303
/* if we accept the proposed network */
305-
if (cutoff >= 0.0 || log(unif_rand()) < cutoff) {
304+
if (cutoff >= 0.0 || log(unif_rand()) < cutoff) {
306305
/* Hold off updating timesteamps until the changes are committed,
307306
which doesn't happen until later. */
308307
for (unsigned int i=0; i < MHp->ntoggles; i++){
@@ -319,17 +318,17 @@ MCMCDynStatus MCMCDyn1Step(ErgmState *s,
319318
sw++; si += i;
320319
sw2 *= sdecay*sdecay; si2 *= sdecay;
321320
sw2++; si2 += i*i;
322-
323-
if(step >= min_MH_interval && !finished) {
321+
322+
if(step >= min_MH_interval && !finished) {
324323
// Now, perform the test:
325324
double mi = (double)si / sw, mi2 = (double)si2 / sw;
326-
325+
327326
double vi = mi2 - mi*mi;
328327
double zi = mi / sqrt(vi * sw2/(sw*sw)); // denom = sqrt(sum(w^2 * v)/sum(w)^2)
329328
double pi = pnorm(zi, 0, 1, FALSE, FALSE); // Pr(Z > zi)
330329

331330
if(verbose>=5) Rprintf("%u: sw=%2.2f sw2=%2.2f d=%d i=%d si=%2.2f si2=%2.2f mi=%2.2f vi=%2.2f ni=%2.2f zi=%2.2f pi=%2.2f\n", step, sw, sw2, kh_size(discord), i, si, si2, mi, vi, (sw*sw)/sw2, zi, pi);
332-
331+
333332
if(pi > MH_pval){
334333
extrasteps = step*MH_interval_add+round(runif(0,1));
335334
finished++;
@@ -340,7 +339,7 @@ MCMCDynStatus MCMCDyn1Step(ErgmState *s,
340339
}
341340

342341
/* Step finished: record changes. */
343-
342+
344343
if(verbose>=4){
345344
if(step>=max_MH_interval ) Rprintf("Convergence not achieved after %u M-H steps.\n",step);
346345
else Rprintf("Convergence achieved after %u M-H steps.\n",step);
@@ -362,14 +361,14 @@ MCMCDynStatus MCMCDyn1Step_advance(ErgmState *s,
362361
int verbose){
363362
StoreDyadMapInt *discord = dur_inf->discord;
364363
int t = dur_inf->time;
365-
364+
366365
Network *nwp = s->nwp;
367366
Model *m = s->m;
368367
MHProposal *MHp = s->MHp;
369368

370369
if(nextdiffedge) {
371370
TailHead dyad;
372-
kh_foreach_key(discord, dyad,{
371+
kh_foreach_key(discord, dyad,{
373372
if(*nextdiffedge<maxchanges){
374373
// and record the toggle.
375374
if(difftime) kv_push(int, *difftime, t);

0 commit comments

Comments
 (0)