[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