[Rsiena-commits] r127 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSienaTest RSienaTest/inst/doc RSienaTest/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 25 18:39:16 CET 2010
Author: ripleyrm
Date: 2010-11-25 18:39:16 +0100 (Thu, 25 Nov 2010)
New Revision: 127
Added:
pkg/RSiena/R/initializeFRAN.r
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/changeLog
pkg/RSiena/inst/doc/s_man400.pdf
pkg/RSiena/man/RSiena-package.Rd
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/changeLog
pkg/RSienaTest/inst/doc/s_man400.pdf
pkg/RSienaTest/man/RSiena-package.Rd
Log:
Forgot initializeFRAN for RSiena and updated scripts for sienaTimeTest in manual
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2010-11-25 17:15:33 UTC (rev 126)
+++ pkg/RSiena/DESCRIPTION 2010-11-25 17:39:16 UTC (rev 127)
@@ -1,7 +1,7 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.126
+Version: 1.0.12.127
Date: 2010-11-25
Author: Various
Depends: R (>= 2.9.0), xtable
Added: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r (rev 0)
+++ pkg/RSiena/R/initializeFRAN.r 2010-11-25 17:39:16 UTC (rev 127)
@@ -0,0 +1,1720 @@
+#/******************************************************************************
+# * 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]
+ }
+ }
+ ## 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)
+ }
+ ## 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]
+ interactions <- effects$effectNumber %in%
+ interactionNos
+ interactionMainEffects <- effects[interactions, ]
+ effects$requested <- effects$include
+ requestedEffects <- effects[effects$include, ]
+
+ effects$include[interactions] <- TRUE
+ effects <- effects[effects$include, ]
+
+ ## split and rejoin both versions before continuing
+ depvarnames <- names(data[[1]]$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"]
+
+ ## 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, "periodNos") <- attr(data, "periodNos")
+ 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
+ ## scores are now available (30.10.10) but still not maxlike.
+ ## check model type: default for symmetric is type 2 (forcing model).
+ syms <- attr(data,"symmetric")
+ z$FinDiffBecauseSymmetric <- FALSE
+ z$modelType <- x$modelType
+ 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 (x$modelType == 1)
+ {
+ z$modelType <- 2
+ }
+ }
+ 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
+ }
+ 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"))
+ if (!all(attr(f,"netnames") %in% effects$name))
+ {
+ stop("Must have at least one effect for each dependent variable")
+ }
+ myeffects <- split(effects, splitFactor)
+ myCompleteEffects <- myeffects
+ ## remove interaction effects and save till later
+ basicEffects <-
+ lapply(myeffects, function(x)
+ {
+ x[!x$shortName %in% c("unspInt", "behUnspInt"), ]
+ }
+ )
+ basicEffectsl <-
+ lapply(myeffects, function(x)
+ {
+ !x$shortName %in% c("unspInt", "behUnspInt")
+ }
+ )
+
+ interactionEffects <-
+ lapply(myeffects, function(x)
+ {
+ x[x$shortName %in% c("unspInt", "behUnspInt"), ]
+ }
+ )
+ interactionEffectsl <-
+ lapply(myeffects, function(x)
+ {
+ x$shortName %in% c("unspInt", "behUnspInt")
+ }
+ )
+ ## store effects objects as we may need to recreate them
+ f$interactionEffects <- interactionEffects
+ f$basicEffects <- basicEffects
+ f$interactionEffectsl <- interactionEffectsl
+ f$basicEffectsl <- basicEffectsl
+ }
+ else
+ {
+ myCompleteEffects <- ff$myCompleteEffects
+ basicEffects <- ff$basicEffects
+ interactionEffects <- ff$interactionEffects
+ basicEffectsl <- ff$basicEffectsl
+ interactionEffectsl <- ff$interactionEffectsl
+ 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 insert in
+ ## effects object in the same rows for later use
+ 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
+ }
+ myCompleteEffects[[i]][basicEffectsl[[i]], ] <- basicEffects[[i]]
+ myCompleteEffects[[i]][interactionEffectsl[[i]],] <-
+ interactionEffects[[i]]
+ ##myeffects[[i]] <- myeffects[[i]][order(myeffects[[i]]$effectNumber),]
+ }
+ ## remove the effects only created as underlying effects
+ ## for interaction effects. first store the original for use next time
+ myeffects <- lapply(myCompleteEffects, 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
+ f$myCompleteEffects <- myCompleteEffects
+ 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$periodNos <- attr(data, "periodNos")
+ 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
+ }
+ else if (ccOption == 3)
+ {
+ ## add the row and column to the missing data
+ mat2[j, ] <- 1
+ mat2[, j] <- 1
+ mat2[j, j] <- 0
+ ## set to missing in raw data for distances later
+ depvar[[i]][j, ] <- NA
+ depvar[[i]][, j] <- NA
+ }
+ }
+ for (j in threes) ## False data is preceded and followed by real
+ {
+ if (ccOption %in% c(1, 2))
+ {
+ ## find missing values for this actor
+ use <- mat2[j, ] > 0
+ ## remove these from mat2, the missing data
+ mat2[j, use] <- 0
+ mat2[use, j] <- 0
+ ## carry forward
+ if (i == 1)
+ {
+ ## 0 any matches from mat1, the real data
+ mat1[j, use] <- 0
+ mat1[use, j] <- 0
+ }
+ else
+ {
+ mat1[j, use] <- networks[[i-1]][j, use]
+ mat1[use, j] <- networks[[i-1]][use, j]
+ }
+ depvar[[i]][j, use] <- 0 ## not missing
+ depvar[[i]][use, j] <- 0
+ depvar[[i]][j, j] <- NA
+ }
+ else if (ccOption == 3)
+ {
+ ## add the row and column to the missing data
+ mat2[j, ] <- 1
+ mat2[, j] <- 1
+ mat2[j, j] <- 0
+ depvar[[i]][j, ] <- NA
+ depvar[[i]][, j] <- NA
+ }
+ }
+ for (j in twos) ## False data is not followed by anything real
+ {
+ if (ccOption == 1)
+ {
+ ## find missing values for this actor
+ use <- mat2[j, ] > 0
+ ## remove these from mat2, the missing data
+ mat2[j, use] <- 0
+ mat2[use, j] <- 0
+ depvar[[i]][j, use] <- 0 ## not missing
+ depvar[[i]][use, j] <- 0
+ depvar[[i]][j, j] <- NA
+ ## carry forward
+ if (i == 1)
+ {
+ ## 0 any matches from mat1, the real data
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 127
More information about the Rsiena-commits
mailing list