[Rsiena-commits] r117 - in pkg: RSiena RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/inst/examples RSienaTest/man RSienaTest/src RSienaTest/src/data RSienaTest/src/model RSienaTest/src/model/ml RSienaTest/src/model/variables
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 20 18:12:15 CEST 2010
Author: ripleyrm
Date: 2010-08-20 18:12:14 +0200 (Fri, 20 Aug 2010)
New Revision: 117
Added:
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/src/siena07internals.cpp
pkg/RSienaTest/src/siena07internals.h
pkg/RSienaTest/src/siena07models.cpp
pkg/RSienaTest/src/siena07models.h
pkg/RSienaTest/src/siena07setup.cpp
pkg/RSienaTest/src/siena07setup.h
pkg/RSienaTest/src/siena07utilities.cpp
pkg/RSienaTest/src/siena07utilities.h
Removed:
pkg/RSienaTest/src/siena07.cpp
Modified:
pkg/RSiena/changeLog
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/R/bayes.r
pkg/RSienaTest/R/maxlikec.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/robmon.r
pkg/RSienaTest/R/siena07.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/RSiena.bib
pkg/RSienaTest/doc/Siena_algorithms4.tex
pkg/RSienaTest/doc/s_man400.tex
pkg/RSienaTest/doc/simstats0c.tex
pkg/RSienaTest/inst/doc/s_man400.pdf
pkg/RSienaTest/inst/examples/algorithms.r
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/src/Makevars
pkg/RSienaTest/src/Makevars.win
pkg/RSienaTest/src/data/NetworkLongitudinalData.cpp
pkg/RSienaTest/src/data/NetworkLongitudinalData.h
pkg/RSienaTest/src/model/EpochSimulation.cpp
pkg/RSienaTest/src/model/EpochSimulation.h
pkg/RSienaTest/src/model/Model.cpp
pkg/RSienaTest/src/model/Model.h
pkg/RSienaTest/src/model/State.cpp
pkg/RSienaTest/src/model/State.h
pkg/RSienaTest/src/model/ml/Chain.cpp
pkg/RSienaTest/src/model/ml/Chain.h
pkg/RSienaTest/src/model/ml/MLSimulation.cpp
pkg/RSienaTest/src/model/ml/MLSimulation.h
pkg/RSienaTest/src/model/ml/NetworkChange.cpp
pkg/RSienaTest/src/model/variables/BehaviorVariable.cpp
pkg/RSienaTest/src/model/variables/BehaviorVariable.h
pkg/RSienaTest/src/model/variables/DependentVariable.cpp
pkg/RSienaTest/src/model/variables/DependentVariable.h
pkg/RSienaTest/src/model/variables/NetworkVariable.cpp
pkg/RSienaTest/src/model/variables/NetworkVariable.h
Log:
ML, symmetric network models, all still work in progress
Modified: pkg/RSiena/changeLog
===================================================================
--- pkg/RSiena/changeLog 2010-08-20 12:12:03 UTC (rev 116)
+++ pkg/RSiena/changeLog 2010-08-20 16:12:14 UTC (rev 117)
@@ -1,3 +1,24 @@
+2010-08-20 R-forge revision 117 RSienaTest only
+
+ * doc/config.dox: added doxygen configuration file
+ * doc/RSiena.bib, s_man400.tex, Siena_algorithms4.tex,
+ simstats0c.tex: documentation updates
+ * inst/examples/algorithms.r: improvements!
+ * R/initialiseFRAN.r, R/simstatsc.r: removed initialise function
+ to separate source file
+ * src/siena07.cpp deleted
+ * src/siena07internals.cpp, src/siena07models.cpp,
+ src/siena07setup.cpp, src/siena07utilities.cpp,
+ src/siena07internals.h, src/siena07models.h,
+ src/siena07setup.h, src/siena07utilities.h: replacements for
+ siena07.cpp
+ * R/bayes.r, R/maxlikec.r, R/phase1.r, R/robmon.r, R/siena07.r,
+ src/many, man/siena07.Rd : ML changes (in progress)
+ * src/model/variables/DependentVariable.cpp,
+ src/model/variables/NetworkVariable.cpp,
+ src/model/variables/DependentVariable.h,
+ src/model/variables/NetworkVariable.h: symmetric networks (in progress)
+
2010-08-20 R-forge revision 116 RSiena only
* R/simstatsc.r: forgotten part of previous change to print.
Modified: pkg/RSienaTest/DESCRIPTION
===================================================================
--- pkg/RSienaTest/DESCRIPTION 2010-08-20 12:12:03 UTC (rev 116)
+++ pkg/RSienaTest/DESCRIPTION 2010-08-20 16:12:14 UTC (rev 117)
@@ -1,7 +1,7 @@
Package: RSienaTest
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.115
+Version: 1.0.12.117
Date: 2010-08-20
Author: Various
Depends: R (>= 2.9.0), xtable
Modified: pkg/RSienaTest/R/bayes.r
===================================================================
--- pkg/RSienaTest/R/bayes.r 2010-08-20 12:12:03 UTC (rev 116)
+++ pkg/RSienaTest/R/bayes.r 2010-08-20 16:12:14 UTC (rev 117)
@@ -24,9 +24,9 @@
z$candidates <- matrix(NA, nrow=nmain * nrunMHBatches,
ncol=sum(!basicRate))
z$acceptances <- rep(NA, nmain * nrunMHBatches)
- z$MHacceptances <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=6)
- z$MHrejections <- matrix(NA, nrow=nmain * nrunMHBatches , ncol=6)
- z$MHproportions <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=6)
+ z$MHacceptances <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=7)
+ z$MHrejections <- matrix(NA, nrow=nmain * nrunMHBatches , ncol=7)
+ z$MHproportions <- matrix(NA, nrow=nmain * nrunMHBatches, ncol=7)
z
}
storeData <- function()
@@ -127,13 +127,13 @@
}
if (z$scaleFactor < tiny)
{
- cat('scalefactor < tiny\n')
+ cat('calefactor < tiny\n')
browser()
}
}
cat('fine tuning took ', iter, ' iterations. Scalefactor:',
z$scaleFactor, '\n')
- z
+ z
}
## initialise
@@ -151,6 +151,7 @@
z$maxlike <- TRUE
model$maxlike <- TRUE
model$FRANname <- "maxlikec"
+ z$print <- FALSE
z$int <- 1
z$int2 <- 1
model$cconditional <- FALSE
@@ -195,7 +196,8 @@
cat('main after ii',ii,numm, '\n')
dev.set(thetaplot)
thetadf <- data.frame(z$lambdas, z$betas)
- acceptsdf <- data.frame(z$MHproportions, z$acceptances)
+ acceptsdf <- data.frame(z$MHproportions[, 1:5],
+ z$acceptances)
lambdaNames <- paste(z$effects$name[basicRate],
z$effects$shortName[basicRate],
z$effects$period[basicRate],
@@ -204,7 +206,8 @@
z$effects$shortName[!basicRate], sep=".")
names(thetadf) <- c(lambdaNames, betaNames)
names(acceptsdf) <- c("InsDiag", "CancDiag", "Permute", "InsPerm",
- "CancPerm", "Missing", "BayesAccepts")
+ "CancPerm", #"Missing",
+ "BayesAccepts")
varnames <- paste(names(thetadf), sep="", collapse= " + ")
varcall <- paste("~ ", varnames, sep="", collapse="")
print(histogram(as.formula(varcall), data=thetadf, scales="free",
@@ -230,7 +233,7 @@
group <- 1
f <- FRANstore()
ans <- .Call("MCMCcycle", PACKAGE=pkgname, f$pData, f$pModel,
- f$pMLSimulation, f$myeffects, as.integer(period),
+ f$myeffects, as.integer(period),
as.integer(group),
z$scaleFactor, nrunMH, nrunMHBatches)
## process the return values
Added: pkg/RSienaTest/R/initializeFRAN.r
===================================================================
--- pkg/RSienaTest/R/initializeFRAN.r (rev 0)
+++ pkg/RSienaTest/R/initializeFRAN.r 2010-08-20 16:12:14 UTC (rev 117)
@@ -0,0 +1,1683 @@
+#/******************************************************************************
+# * SIENA: Simulation Investigation for Empirical Network Analysis
+# *
+# * Web: http://www.stats.ox.ac.uk/~snidjers/siena
+# *
+# * File: initializeFRAN.r
+# *
+# * Description: This module contains the code for setting up the data in C++
+# * and for ML making the initial chain.
+# *****************************************************************************/
+##@initializeFRAN siena07 reformat data and send to C. get targets.
+initializeFRAN <- function(z, x, data, effects, prevAns, initC, profileData,
+ returnDeps)
+{
+ ## fix up the interface so can call from outside robmon framework
+ if (is.null(z$FinDiff.method))
+ {
+ z$FinDiff.method <- FALSE
+ }
+ if (is.null(z$int))
+ {
+ z$int <- 1
+ }
+ if (is.null(z$int2))
+ {
+ z$int2 <- 1
+ }
+
+ if (!initC) ## ie first time round
+ {
+ if (!inherits(data,'siena'))
+ {
+ stop('not valid siena data object')
+ }
+ ## check the effects object
+ defaultEffects <- getEffects(data)
+ if (is.null(effects))
+ {
+ effects <- defaultEffects
+ }
+ else
+ {
+ ## todo check that the effects match the data dependent variables
+ userlist <- apply(effects[effects$include,], 1, function(x)
+ paste(x[c("name", "effectName",
+ "type", "groupName")],
+ collapse="|"))
+ deflist <- apply(defaultEffects, 1, function(x)
+ paste(x[c("name", "effectName",
+ "type", "groupName")],
+ collapse="|"))
+ if (!all(userlist %in% deflist))
+ {
+ bad <- which(!(userlist %in% deflist))
+ print(userlist[bad])
+ stop("invalid effect requested: see above ")
+ }
+ }
+ if (!inherits(effects, 'data.frame'))
+ {
+ stop('effects is not a data.frame')
+ }
+ if (x$useStdInits)
+ {
+ if (any(effects$effectName != defaultEffects$effectName))
+ {
+ stop('Cannot use standard initialisation with a ',
+ 'different effect list')
+ }
+ effects$initialValue <- defaultEffects$initialValue
+ }
+ else
+ {
+ if (!is.null(prevAns) && inherits(prevAns, "sienaFit"))
+ {
+ prevEffects <- prevAns$requestedEffects
+ prevEffects$initialValue <- prevAns$theta
+ if (prevAns$cconditional)
+ {
+ condEffects <- attr(prevAns$f, "condEffects")
+ condEffects$initialValue <- prevAns$rate
+ prevEffects <- rbind(prevEffects, condEffects)
+ }
+ oldlist <- apply(prevEffects, 1, function(x)
+ paste(x[c("name", "shortName",
+ "type", "groupName",
+ "interaction1", "interaction2",
+ "period")],
+ collapse="|"))
+ efflist <- apply(effects, 1, function(x)
+ paste(x[c("name", "shortName",
+ "type", "groupName",
+ "interaction1", "interaction2",
+ "period")],
+ collapse="|"))
+ use <- efflist %in% oldlist
+ effects$initialValue[use] <-
+ prevEffects$initialValue[match(efflist, oldlist)][use]
+ }
+ }
+ ## add any effects needed for time dummies
+ tmp <- sienaTimeFix(effects, data)
+ data <- tmp$data
+ effects <- tmp$effects
+ ## find any effects not included which are needed for interactions
+ interactionNos <- unique(c(effects$effect1, effects$effect2,
+ effects$effect3))
+ interactionNos <- interactionNos[interactionNos > 0]
+ interactionMainEffects <- effects[interactionNos, ]
+ effects$requested <- effects$include
+ requestedEffects <- effects[effects$include, ]
+
+ effects$include[interactionNos] <- TRUE
+ effects <- effects[effects$include,]
+
+ ## split and rejoin both versions before continuing
+ if (inherits(data, "sienaGroup"))
+ depvarnames <- names(data[[1]]$depvars)
+ else
+ depvarnames <- names(data$depvars)
+
+ effects1 <- split(requestedEffects, requestedEffects$name)
+ effects1order <- match(depvarnames, names(effects1))
+ requestedEffects <- do.call(rbind, effects1[effects1order])
+ row.names(requestedEffects) <- 1:nrow(requestedEffects)
+
+ effects1 <- split(effects, effects$name)
+ effects1order <- match(depvarnames, names(effects1))
+ effects <- do.call(rbind, effects1[effects1order])
+ row.names(effects) <- 1:nrow(effects)
+
+ ## now set up z provisionally
+ z$theta <- requestedEffects$initialValue
+ z$fixed <- requestedEffects$fix
+ z$test <- requestedEffects$test
+ z$pp <- length(z$test)
+ z$posj <- rep(FALSE, z$pp)
+ z$posj[requestedEffects$basicRate] <- TRUE
+ z$BasicRateFunction <- z$posj
+
+ ## sort out names of user specified interaction effects
+ effects <- fixUpEffectNames(effects)
+ ## copy interaction names to the requested effects
+ requestedEffects$effectName <- effects[effects$requested,
+ "effectName"]
+ requestedEffects$functionName <- effects[effects$requested,
+ "functionName"]
+
+ ## get data object into group format to save coping with two
+ ## different formats
+ if (inherits(data, 'sienaGroup'))
+ {
+ nGroup <- length(data)
+ }
+ else
+ {
+ nGroup <- 1
+ data <- sienaGroupCreate(list(data), singleOK=TRUE)
+ }
+ ## if not specified whether conditional or nor, set to conditional
+ ## iff there is only one dependent variable (therefore number 1)
+ ## and not maxlike
+ if (is.na(x$cconditional))
+ {
+ x$cconditional <- !x$maxlike && (length(depvarnames) == 1)
+ if (x$cconditional)
+ {
+ x$condvarno <- 1
+ }
+ }
+ types <- sapply(data[[1]]$depvars, function(x) attr(x, 'type'))
+ nets <- sum(types != "behavior")
+ ## now check if conditional estimation is OK and copy to z if so
+ z$cconditional <- FALSE
+ if (x$cconditional)
+ {
+ if (x$maxlike)
+ {
+ stop("Conditional estimation is not possible with",
+ "maximum likelihood method")
+ }
+ ## if (nets == 1) not sure if this is necessary
+ ## {
+ z$cconditional <- TRUE
+ ## find the conditioning variable
+ observations <- attr(data, "observations")
+ ## this is actual number of waves to process
+ if (x$condname != "")
+ {
+ z$condvarno <- match(x$condname, attr(data, "netnames"))
+ z$condname <- x$condname
+ }
+ else
+ {
+ z$condvarno <- x$condvarno
+ z$condname <- attr(data, "netnames")[x$condvarno]
+ }
+ z$condtype <- attr(data, "types")[z$condvarno]
+ if (z$condtype == "oneMode")
+ z$symmetric <- attr(data, "symmetric")[[z$condvarno]]
+ else
+ z$symmetric <- FALSE
+ ## find the positions of basic rate effects for this network
+ z$condvar <-
+ (1:nrow(requestedEffects))[requestedEffects$name==
+ z$condname][1:observations]
+ z$theta<- z$theta[-z$condvar]
+ z$fixed<- z$fixed[-z$condvar]
+ z$test<- z$test[-z$condvar]
+ z$pp<- z$pp-length(z$condvar)
+ z$scale<- z$scale[-z$condvar]
+ z$BasicRateFunction <- z$posj[-z$condvar]
+ z$posj <- z$posj[-z$condvar]
+ z$theta[z$posj] <-
+ z$theta[z$posj] / requestedEffects$initialValue[z$condvar]
+ }
+
+ ## unpack data and put onto f anything we may need next time round.
+ f <- lapply(data, function(x) unpackData(x))
+
+ attr(f, "netnames") <- attr(data, "netnames")
+ attr(f, "symmetric") <- attr(data, "symmetric")
+ attr(f, "allUpOnly") <- attr(data, "allUpOnly")
+ attr(f, "allDownOnly") <- attr(data, "allDownOnly")
+ attr(f, "allHigher") <- attr(data, "allHigher")
+ attr(f, "allDisjoint") <- attr(data, "allDisjoint")
+ attr(f, "allAtLeastOne") <- attr(data, "allAtLeastOne")
+ attr(f, "anyUpOnly") <- attr(data, "anyUpOnly")
+ attr(f, "anyDownOnly") <- attr(data, "anyDownOnly")
+ attr(f, "anyHigher") <- attr(data, "anyHigher")
+ attr(f, "anyDisjoint") <- attr(data, "anyDisjoint")
+ attr(f, "anyAtLeastOne") <- attr(data, "anyAtLeastOne")
+ attr(f, "types") <- attr(data, "types")
+ attr(f, "observations") <- attr(data, "observations")
+ attr(f, "compositionChange") <- attr(data, "compositionChange")
+ attr(f, "exooptions") <- attr(data, "exooptions")
+ attr(f, "groupPeriods") <- attr(data, "groupPeriods")
+ attr(f, "numberNonMissingNetwork") <-
+ attr(data, "numberNonMissingNetwork")
+ attr(f, "numberMissingNetwork") <- attr(data, "numberMissingNetwork")
+ attr(f, "numberNonMissingBehavior") <-
+ attr(data, "numberNonMissingBehavior")
+ attr(f, "numberMissingBehavior") <- attr(data, "numberMissingBehavior")
+
+ # attr(f, "totalMissings") <- attr(data, "totalMissings")
+
+ if (x$maxlike && x$FinDiff.method)
+ {
+ stop("Finite difference method for derivatives not available",
+ "with Maximum likelihood method")
+ }
+ ## if any networks symmetric must use finite differences and not maxlike
+ syms <- attr(data,"symmetric")
+ # z$FinDiffBecauseSymmetric <- FALSE
+ # if (any(!is.na(syms) & syms))
+ # {
+ # z$FinDiff.method <- TRUE
+ # z$FinDiffBecauseSymmetric <- TRUE
+ # if (x$maxlike)
+ # {
+ # stop("Maximum likelihood method not implemented",
+ # "for symmetric networks")
+ # }
+ # }
+ if (z$cconditional)
+ {
+ attr(f, "change") <-
+ sapply(f, function(xx)attr(xx$depvars[[z$condname]],
+ 'distance'))
+ attr(f,"condEffects") <- requestedEffects[z$condvar,]
+ effcondvar <-
+ (1:nrow(effects))[effects$name==
+ z$condname][1:observations]
+ effects <- effects[-effcondvar, ]
+ requestedEffects <- requestedEffects[-z$condvar,]
+ }
+ ## use previous dfra only if everything matches
+ if (!is.null(prevAns) && inherits(prevAns, "sienaFit"))
+ {
+ if ((nrow(prevAns$dfra) == nrow(requestedEffects)) &&
+ all(rownames(prevAns$dfra) == requestedEffects$shortName)
+ && !is.null(prevAns$sf))
+ {
+ z$haveDfra <- TRUE
+ z$dfra <- prevAns$dfra
+ z$dinv <- prevAns$dinv
+ z$sf <- prevAns$sf
+ }
+ }
+ z$effects <- effects
+ z$requestedEffects <- requestedEffects
+ }
+ else ## initC, i.e just send already set up data into new processes
+ {
+ f <- FRANstore()
+ ## Would like f to be just the data objects plus the attributes
+ ## but need the effects later. Also a few other things,
+ ## which probably could be attributes but are not!
+ ## They will be automatically removed: if needed they must be readded
+ ff <- f
+ nGroup <- f$nGroup
+ f[(nGroup + 1): length(f)] <- NULL
+ }
+ ##browser()
+ pData <- .Call('setupData', PACKAGE=pkgname,
+ lapply(f, function(x)(as.integer(x$observations))),
+ lapply(f, function(x)(x$nodeSets)))
+ ans <- .Call('OneMode', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$nets))
+ ans <- .Call('Bipartite', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$bipartites))
+ ans <- .Call('Behavior', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$behavs))
+ ans <-.Call('ConstantCovariates', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$cCovars))
+ ans <-.Call('ChangingCovariates', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$vCovars))
+ ans <-.Call('DyadicCovariates', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$dycCovars))
+ ans <-.Call('ChangingDyadicCovariates', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$dyvCovars))
+ ans <-.Call('ExogEvent', PACKAGE=pkgname,
+ pData, lapply(f, function(x)x$exog))
+
+ ## split the names of the constraints
+ higher <- attr(f, "allHigher")
+ disjoint <- attr(f, "allDisjoint")
+ atLeastOne <- attr(f, "allAtLeastOne")
+ froms <- sapply(strsplit(names(higher), ","), function(x)x[1])
+ tos <- sapply(strsplit(names(higher), ","), function(x)x[2])
+ ans <- .Call("Constraints", PACKAGE=pkgname,
+ pData, froms[higher], tos[higher],
+ froms[disjoint], tos[disjoint],
+ froms[atLeastOne], tos[atLeastOne])
+
+ ##store the address
+ f$pData <- pData
+ ## register a finalizer
+ ans <- reg.finalizer(f$pData, clearData, onexit = FALSE)
+
+ if (!initC)
+ {
+ storage.mode(effects$parm) <- 'integer'
+ storage.mode(effects$group) <- 'integer'
+ storage.mode(effects$period) <- 'integer'
+ effects$effectPtr <- rep(NA, nrow(effects))
+ splitFactor <- factor(effects$name, levels=attr(f, "netnames"))
+ myeffects <- split(effects, splitFactor)
+ ## remove interaction effects and save till later
+ basicEffects <-
+ lapply(myeffects, function(x)
+ {
+ x[!x$shortName %in% c("unspInt", "behUnspInt"), ]
+ }
+ )
+ interactionEffects <-
+ lapply(myeffects, function(x)
+ {
+ x[x$shortName %in% c("unspInt", "behUnspInt"), ]
+ }
+ )
+ ## store effects objects as we may need to recreate them
+ f$interactionEffects <- interactionEffects
+ f$basicEffects <- basicEffects
+ }
+ else
+ {
+ myeffects <- ff$myeffects
+ basicEffects <- ff$basicEffects
+ interactionEffects <- ff$interactionEffects
+ types <- ff$types
+ }
+ ans <- .Call('effects', PACKAGE=pkgname, pData, basicEffects)
+ pModel <- ans[[1]][[1]]
+ for (i in seq(along=(ans[[2]]))) ## ans[[2]] is a list of lists of
+ ## pointers to effects. Each list corresponds to one
+ ## dependent variable
+ {
+ effectPtr <- ans[[2]][[i]]
+ basicEffects[[i]]$effectPtr <- effectPtr
+ interactionEffects[[i]]$effect1 <-
+ basicEffects[[i]]$effectPtr[match(interactionEffects[[i]]$effect1,
+ basicEffects[[i]]$effectNumber)]
+ interactionEffects[[i]]$effect2 <-
+ basicEffects[[i]]$effectPtr[match(interactionEffects[[i]]$effect2,
+ basicEffects[[i]]$effectNumber)]
+ interactionEffects[[i]]$effect3 <-
+ basicEffects[[i]]$effectPtr[match(interactionEffects[[i]]$effect3,
+ basicEffects[[i]]$effectNumber)]
+ }
+ ans <- .Call('interactionEffects', PACKAGE=pkgname,
+ pData, pModel, interactionEffects)
+ ## copy these pointers to the interaction effects and then rejoin
+ for (i in 1:length(ans[[1]])) ## ans is a list of lists of
+ ## pointers to effects. Each list corresponds to one
+ ## dependent variable
+ {
+ if (nrow(interactionEffects[[i]]) > 0)
+ {
+ effectPtr <- ans[[1]][[i]]
+ interactionEffects[[i]]$effectPtr <- effectPtr
+ }
+ myeffects[[i]] <- rbind(basicEffects[[i]], interactionEffects[[i]])
+ }
+ ## remove the effects only created as underlying effects
+ ## for interaction effects
+ myeffects <- lapply(myeffects, function(x)
+ {
+ x[x$requested, ]
+ }
+ )
+ if (!initC)
+ {
+ ans <- .Call("getTargets", PACKAGE=pkgname, pData, pModel, myeffects)
+ if (!x$maxlike)
+ {
+ z$targets <- rowSums(ans)
+ z$targets2 <- ans
+ }
+ else
+ {
+ z$targets <- rep(0, z$pp)
+ z$targets2 <- 0
+ z$maxlikeTargets <- rowSums(ans)
+ z$maxlikeTargets2 <- ans
+ z$mult <- 1
+ z$nrunMH <- z$mult * sum(z$maxlikeTargets[z$effects$basicRate])
+ }
+ }
+
+ ##store address of model
+ f$pModel <- pModel
+ ans <- reg.finalizer(f$pModel, clearModel, onexit = FALSE)
+ if (x$MaxDegree == 0 || is.null(x$MaxDegree))
+ {
+ MAXDEGREE <- NULL
+ }
+ else
+ {
+ MAXDEGREE <- as.integer(x$MaxDegree)
+ storage.mode(MAXDEGREE) <- "integer"
+ }
+ if (z$cconditional)
+ {
+ CONDVAR <- z$condname
+ CONDTARGET <- attr(f, "change")
+ ## cat(CONDTARGET, '\n')
+ }
+ else
+ {
+ CONDVAR <- NULL
+ CONDTARGET <- NULL
+ }
+
+ ans <- .Call("setupModelOptions", PACKAGE=pkgname,
+ pData, pModel, MAXDEGREE, CONDVAR, CONDTARGET,
+ profileData, z$parallelTesting, x$MODELTYPE)
+ if (x$maxlike)
+ {
+ ## set up chains and do initial steps
+ simpleRates <- TRUE
+
+ types <- attr(f, "types")
+ nbrNonMissNet <- attr(f, "numberNonMissingNetwork")
+ nbrMissNet <- attr(f, "numberMissingNetwork")
+ nbrNonMissBeh <- attr(f, "numberNonMissingBehavior")
+ nbrMissBeh <- attr(f, "numberMissingBehavior")
+
+ z$prmin <- nbrMissNet/ (nbrMissNet + nbrNonMissNet)
+ if (sum(nbrMissBeh + nbrNonMissBeh) > 0)
+ {
+ z$prmib <- nbrMissBeh/ (nbrMissBeh + nbrNonMissBeh)
+ }
+ else
+ {
+ z$prmib <- rep(0, length(nbrMissBeh))
+ }
+ ## cat (z$prmin, z$prmib, '\n')
+ z$probs <- c(x$pridg, x$prcdg, x$prper, x$pripr, x$prdpr, x$prirms,
+ x$prdrms)
+ cat(z$probs,'\n')
+ ans <- .Call("mlMakeChains", PACKAGE=pkgname, pData, pModel,
+ simpleRates, z$probs, z$prmin, z$prmib,
+ x$minimumPermutationLength,
+ x$maximumPermutationLength,
+ x$initialPermutationLength)
+
+ f$minimalChain <- ans[[1]]
+ f$chain <- ans[[2]]
+ # print(nrow(ans[[1]][[1]]))
+ # print(nrow(ans[[2]][[1]]))
+ # browser()
+ ##store address of simulation object
+ # f$pMLSimulation <- ans[[1]][[1]]
+ #ans1 <- reg.finalizer(f$pMLSimulation, clearMLSimulation,
+ # onexit = FALSE)
+ }
+ f$myeffects <- myeffects
+ if (!initC)
+ {
+ if (is.null(z$print) || z$print)
+ {
+ DataReport(z, x, f)
+ }
+ f$randomseed2 <- z$randomseed2
+ }
+ else
+ {
+ f$randomseed2 <- ff$randomseed2
+ }
+ f$observations <- attr(f, "observations") + 1
+ f$depNames <- names(f[[1]]$depvars)
+ f$groupNames <- names(f)[1:nGroup]
+ f$nGroup <- nGroup
+ f$types <- types
+ if (!initC)
+ {
+ z$f <- f
+ z <- initForAlgorithms(z)
+ z$f[1:nGroup] <- NULL
+ }
+ if (initC || (z$int == 1 && z$int2 == 1 &&
+ (is.null(z$nbrNodes) || z$nbrNodes == 1)))
+ {
+
+ f[1:nGroup] <- NULL
+ }
+ FRANstore(f) ## store f in FRANstore
+ z$returnDeps <- returnDeps
+ z$returnDepsStored <- returnDeps
+ z$observations <- f$observations
+ z
+}
+##@createEdgeLists siena07 Reformat data for C++
+createEdgeLists<- function(mat, matorig)
+{
+ ## mat1 is basic values, with missings and structurals replaced
+ tmp <- lapply(1 : nrow(mat), function(x, y)
+ {
+ mymat <- matrix(0, nrow = sum(y[x, ] > 0), ncol = 3)
+ mymat[, 1] <- x
+ mymat[, 2] <- which(y[x, ] != 0)
+ mymat[, 3] <- y[x, mymat[, 2]]
+ mymat
+ }, y = mat)
+ mat1 <- do.call(rbind, tmp)
+ ## mat2 reverts to matorig to get the missing values
+ tmp <- lapply(1 : nrow(matorig), function(x, y)
+ {
+ mymat <- matrix(0, nrow = sum(is.na(y[x, ])), ncol = 3)
+ mymat[, 1] <- x
+ mymat[, 2] <- which(is.na(y[x, ]))
+ mymat[, 3] <- 1
+ mymat
+ }, y = matorig)
+ mat2 <- do.call(rbind, tmp)
+ ## remove the diagonal
+ mat2 <- mat2[mat2[, 1] != mat2[, 2], , drop=FALSE]
+ ## mat3 structurals
+ struct <- mat1[,3] %in% c(10, 11)
+ mat1[struct, 3] <- mat1[struct,3] - 10
+ mat3 <- mat1[struct, , drop=FALSE]
+ mat3[, 3] <- 1
+ mat1 <- mat1[!mat1[,3] == 0, , drop=FALSE] ##remove any zeros just created
+ ##fix up storage mode to be integer
+ storage.mode(mat1) <- 'integer'
+ storage.mode(mat2) <- 'integer'
+ storage.mode(mat3) <- 'integer'
+ ## add attribute of size
+ nodeSets <- attr(matorig, "nodeSet")
+ if (length(nodeSets) > 1) ## bipartite
+ {
+ attr(mat1,'nActors') <- c(nrow(mat), ncol(mat))
+ attr(mat2,'nActors') <- c(nrow(mat), ncol(mat))
+ attr(mat3,'nActors') <- c(nrow(mat), ncol(mat))
+ }
+ else
+ {
+ attr(mat1,'nActors') <- nrow(mat)
+ attr(mat2,'nActors') <- nrow(mat)
+ attr(mat3,'nActors') <- nrow(mat)
+ }
+
+ list(mat1 = t(mat1), mat2 = t(mat2), mat3 = t(mat3))
+}
+##@createCovarEdgeLists siena07 Reformat data for C++
+createCovarEdgeList<- function(mat, matorig)
+{
+ tmp <- lapply(1 : nrow(mat), function(x, y)
+ {
+ mymat <- matrix(0, nrow = sum(y[x, ] != 0), ncol = 3)
+ mymat[, 1] <- x
+ mymat[, 2] <- which(y[x, ] != 0)
+ mymat[, 3] <- y[x, mymat[, 2]]
+ mymat
+ }, y = mat)
+ mat1 <- do.call(rbind, tmp)
+ ##mat2 reverts to matorig to get the missing values
+ tmp <- lapply(1 : nrow(matorig), function(x, y)
+ {
+ mymat <- matrix(0, nrow = sum(is.na(y[x, ])), ncol = 3)
+ mymat[, 1] <- x
+ mymat[, 2] <- which(is.na(y[x, ]))
+ mymat[, 3] <- 1
+ mymat
+ }, y = matorig)
+ mat2 <- do.call(rbind, tmp)
+ ## add attribute of size
+ attr(mat1,'nActors1') <- nrow(mat)
+ attr(mat1,'nActors2') <- ncol(mat)
+ list(mat1=t(mat1), mat2=t(mat2))
+}
+##@unpackOneMode siena07 Reformat data for C++
+unpackOneMode <- function(depvar, observations, compositionChange)
+{
+ edgeLists <- vector('list', observations)
+ networks <- vector('list', observations)
+ actorSet <- attr(depvar, "nodeSet")
+ compActorSets <- sapply(compositionChange, function(x)attr(x, "nodeSet"))
+ thisComp <- match(actorSet, compActorSets)
+ compChange <- !is.na(thisComp)
+ if (compChange)
+ {
+ action <- attr(compositionChange[[thisComp]], "action")
+ ccOption <- attr(compositionChange[[thisComp]], "ccOption")
+ }
+ else
+ {
+ ccOption <- 0
+ action <- matrix(0, nrow=attr(depvar, "netdims")[1], ncol=observations)
+ }
+ ## sort out composition change
+ ## convertToStructuralZeros()?
+ sparse <- attr(depvar, 'sparse')
+ if (sparse)
+ {
+ ## require(Matrix)
+ ## have a list of sparse matrices in triplet format
+ ## with missings and structurals embedded and 0 based indices!
+ netmiss <- vector("list", observations)
+ for (i in 1:observations)
+ {
+ ## extract this matrix
+ networks[[i]] <- depvar[[i]]
+ nActors <- nrow(depvar[[i]])
+ ## stop if any duplicates
+ netmat <- cbind(networks[[i]]@i+1, networks[[i]]@j+1,
+ networks[[i]]@x)
+ if (any(duplicated(netmat[, 1:2])))
+ {
+ stop("duplicate entries in sparse matrix")
+ }
+ ## extract missing entries
+ netmiss[[i]] <- netmat[is.na(netmat[,3]), , drop = FALSE]
+ netmiss[[i]] <-
+ netmiss[[i]][netmiss[[i]][, 1] != netmiss[[i]][, 2], ,
+ drop=FALSE]
+ ## carry forward missing values if any
+ if (i == 1)
+ {
+ netmat <- netmat[!is.na(netmat[,3]), ]
+ networks[[i]] <- spMatrix(nActors, nActors, netmat[, 1],
+ netmat[, 2], netmat[,3])
+ }
+ else
+ {
+ netmiss1 <- netmiss[[i]][, 1:2]
+ storage.mode(netmiss1) <- 'integer'
+ networks[[i]][netmiss1[, 1:2]] <-
+ networks[[i-1]][netmiss1[, 1:2]]
+ }
+ }
+ for (i in 1:observations)
+ {
+ mat1 <- networks[[i]]
+ ## drop the diagonal, if present
+ diag(mat1) <- 0
+ mat1 <- cbind(mat1 at i + 1, mat1 at j + 1, mat1 at x)
+ ##missing edgelist
+ mat2 <- netmiss[[i]]
+ mat2[, 3] <- 1
+ ## rows of mat1 with structural values
+ struct <- mat1[, 3] %in% c(10, 11)
+ ## reset real data
+ mat1[struct, 3] <- mat1[struct, 3] - 10
+ ## copy reset data to structural edgelist
+ mat3 <- mat1[struct, , drop = FALSE]
+ mat3[, 3] <- 1
+ ## now remove the zeros from reset data
+ mat1 <- mat1[!mat1[, 3] == 0, ]
+ ## do comp change
+ if (compChange)
+ {
+ ## revert to sparse matrices temporarily
+ mat1 <- spMatrix(nrow=nActors, ncol=nActors, i = mat1[, 1],
+ j=mat1[, 2], x=mat1[, 3])
+ mat2 <- spMatrix(nrow=nActors, ncol=nActors, i = mat2[, 1],
+ j=mat2[, 2], x=mat2[, 3])
+ mat3 <- spMatrix(nrow=nActors, ncol=nActors, i = mat3[, 1],
+ j=mat3[, 2], x=mat3[, 3])
+ ones <- which(action[, i] == 1)
+ twos <- which(action[, i] == 2)
+ threes <- which(action[, i] == 3)
+ for (j in ones) ## False data is not preceded by anything real
+ {
+ if (ccOption %in% c(1, 2))
+ {
+ ## find missing values for this actor
+ use <- mat2[j, ] > 0
+ ## remove from real data (i.e. zero)
+ mat1[j, use] <- 0
+ mat1[use, j] <- 0
+ ## remove from missing data
+ mat2[j, use] <- 0
+ mat2[use, j] <- 0
+ ## remove from raw data for distances later
+ depvar[[i]][j, use] <- 0 ## zero
+ depvar[[i]][use, j] <- 0
+ depvar[[i]][j, j] <- NA
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 117
More information about the Rsiena-commits
mailing list