[Rsiena-commits] r77 - in pkg/RSiena: . R data inst/examples man tests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Apr 12 18:38:29 CEST 2010


Author: ripleyrm
Date: 2010-04-12 18:38:28 +0200 (Mon, 12 Apr 2010)
New Revision: 77

Modified:
   pkg/RSiena/R/globals.r
   pkg/RSiena/R/phase1.r
   pkg/RSiena/R/phase2.r
   pkg/RSiena/R/phase3.r
   pkg/RSiena/R/printDataReport.r
   pkg/RSiena/R/printInitialDescription.r
   pkg/RSiena/R/robmon.r
   pkg/RSiena/R/siena01.r
   pkg/RSiena/R/siena07.r
   pkg/RSiena/R/sienaDataCreate.r
   pkg/RSiena/R/sienaDataCreateFromSession.r
   pkg/RSiena/R/sienaModelCreate.r
   pkg/RSiena/R/sienaprint.r
   pkg/RSiena/R/simstatsc.r
   pkg/RSiena/changeLog
   pkg/RSiena/data/allEffects.csv
   pkg/RSiena/inst/examples/EIESATT.dat
   pkg/RSiena/inst/examples/FR0.DAT
   pkg/RSiena/inst/examples/FR1.DAT
   pkg/RSiena/inst/examples/FR2.DAT
   pkg/RSiena/inst/examples/FR3.DAT
   pkg/RSiena/inst/examples/FR4.DAT
   pkg/RSiena/inst/examples/FR5.DAT
   pkg/RSiena/inst/examples/FR6.DAT
   pkg/RSiena/inst/examples/VARS.DAT
   pkg/RSiena/inst/examples/VRND32T0.DAT
   pkg/RSiena/inst/examples/VRND32T1.DAT
   pkg/RSiena/inst/examples/VRND32T2.DAT
   pkg/RSiena/inst/examples/VRND32T5.DAT
   pkg/RSiena/inst/examples/cbc3.dat
   pkg/RSiena/inst/examples/klas12b-advice.dat
   pkg/RSiena/inst/examples/klas12b-alcohol.dat
   pkg/RSiena/inst/examples/klas12b-delinquency.dat
   pkg/RSiena/inst/examples/klas12b-demographics.dat
   pkg/RSiena/inst/examples/klas12b-net-1.dat
   pkg/RSiena/inst/examples/klas12b-net-2.dat
   pkg/RSiena/inst/examples/klas12b-net-3.dat
   pkg/RSiena/inst/examples/klas12b-net-3m.dat
   pkg/RSiena/inst/examples/klas12b-net-4.dat
   pkg/RSiena/inst/examples/klas12b-net-4m.dat
   pkg/RSiena/inst/examples/klas12b-present.dat
   pkg/RSiena/inst/examples/klas12b-primary.dat
   pkg/RSiena/inst/examples/klas12b_netbeh.csv
   pkg/RSiena/inst/examples/klas12b_netbeha.csv
   pkg/RSiena/man/coDyadCovar.Rd
   pkg/RSiena/man/includeEffects.Rd
   pkg/RSiena/man/siena07.Rd
   pkg/RSiena/man/sienaModelCreate.Rd
   pkg/RSiena/tests/parallel.Rout.save
Log:
bug fixes., missing trailing new lines.

Modified: pkg/RSiena/R/globals.r
===================================================================
--- pkg/RSiena/R/globals.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/globals.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -139,9 +139,13 @@
         Report(c("@", level, "\n", text, "\n"), hdest=dest, sep="", fill=fill)
         Report(rep(ch, sum(nchar(text))), hdest=dest, sep="", fill=fill)
         if (level < 3)
+        {
             Report("\n\n", hdest = dest)
+        }
         else
+        {
             Report("\n", hdest = dest)
+        }
     }
 }
 
@@ -149,10 +153,14 @@
 PrtOutMat<- function(mat, dest)
 {
     if (missing(dest))
-        Report(format(t(mat)), sep=c(rep.int(" ", ncol(mat) - 1), "\n"))
+    {
+        Report(format(t(mat), scientific=FALSE),
+               sep=c(rep.int(" ", ncol(mat) - 1), "\n"))
+    }
     else
     {
-        Report(format(t(mat)), sep=c(rep.int(" ", ncol(mat) - 1), "\n"),
+        Report(format(t(mat), scientific=FALSE),
+               sep=c(rep.int(" ", ncol(mat) - 1), "\n"),
                hdest=deparse(substitute(dest)))
         Report("\n", hdest=deparse(substitute(dest)))
     }

Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/phase1.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -46,21 +46,15 @@
     z$sf2 <- array(0, dim=c(z$n1, f$observations - 1, z$pp))
     z$ssc <- array(0, dim=c(z$n1, f$observations - 1, z$pp))
     z$sdf <- array(0, dim=c(z$n1, z$pp, z$pp))
+    z$accepts <- matrix(0, nrow=z$n1, ncol=6)
+    z$rejects <- matrix(0, nrow=z$n1, ncol=6)
     z$npos <- rep(0, z$pp)
     z$writefreq <- 10
     z$DerivativeProblem <- FALSE
     z$Deriv <- !z$FinDiff.method ## can do both in phase 3 but not here!
-    zsmall <- NULL
-    zsmall$theta <- z$theta
-    zsmall$Deriv <- z$Deriv
-    zsmall$Phase <- z$Phase
-    zsmall$nit <- z$nit
-    zsmall$FinDiff.method <- z$FinDiff.method
-    zsmall$int2 <- z$int2
-    zsmall$cl <- z$cl
-    xsmall<- NULL
-    zsmall$cconditional <- z$cconditional
-    zsmall$condvar <- z$condvar
+    xsmall <- NULL
+    zsmall <- makeZsmall(z)
+
     nits <- seq(1, firstNit, int)
     if (any(nits >= 6))
     {
@@ -188,6 +182,7 @@
         z$phase1Its <- z$phase1Its + int
       #  browser()
     }
+   ## browser()
     if (int == 1)
     {
         fra <- colSums(zz$fra)
@@ -208,15 +203,17 @@
         }
 
     }
-    if (z$FinDiff.method)
+    if (x$maxlike)
     {
+        z$sdf[z$nit, , ] <- zz$dff
+        z$accepts[z$nit, ] <- zz$accepts
+        z$rejects[z$nit, ] <- zz$rejects
+    }
+    else if (z$FinDiff.method)
+    {
         z <- FiniteDifferences(z, x, fra + z$targets, ...)
         z$sdf[z$nit:(z$nit + (z$int - 1)), , ] <- z$sdf0
     }
-    else if (x$maxlike)
-    {
-        z$sdf[z$nit, , ] <- zz$dff
-    }
     else
     {
         if (int==1)
@@ -271,17 +268,8 @@
 phase1.2 <- function(z, x, ...)
 {
     ##finish phase 1 iterations and do end-of-phase processing
-    zsmall <- NULL
-    zsmall$theta <- z$theta
-    zsmall$Deriv <- z$Deriv
-    zsmall$Phase <- z$Phase
-    zsmall$nit <- z$nit
-    zsmall$int2 <- z$int2
-    zsmall$cl <- z$cl
-    xsmall<- NULL
-    zsmall$cconditional <- z$cconditional
-    zsmall$condvar <- z$condvar
-    zsmall$FinDiff.method <- z$FinDiff.method
+    zsmall <- makeZsmall(z)
+    xsmall <- NULL
     int <- z$int
     if (z$n1 > z$phase1Its)
     {
@@ -485,12 +473,13 @@
 {
     int <- z$int
     fras <- array(0, dim = c(int, z$pp, z$pp))
-    xsmall<- NULL
+    xsmall <- NULL
  ##browser()
     for (i in 1 : z$pp)
     {
-        zdummy <- z[c('theta', 'Deriv', 'cconditional', 'FinDiff.method',
-                      'int2', 'cl')]
+       # zdummy <- z[c('theta', 'Deriv', 'cconditional', 'FinDiff.method',
+       #               'int2', 'cl')]
+        zdummy <- makeZsmall(z)
         if (!z$fixed[i])
         {
             zdummy$theta[i] <- z$theta[i] + z$epsilon[i]
@@ -522,7 +511,7 @@
                 fras[j, i, ] <- colSums(zz[[j]]$fra) - fra
         }
     }
-                                        # browser()
+                                        ##browser()
     if (z$Phase == 1 && z$nit <= 10)
     {
         for (ii in 1: min(10 - z$nit + 1, int))
@@ -537,7 +526,7 @@
         sdf[i, , ] <- fras[i, , ] / rep(z$epsilon, z$pp)
     }
     z$sdf0 <- sdf
-                                        # browser()
+                                        ## browser()
     z
 }
 ##@derivativeFromScoresAndDeviations siena07 create dfra from scores and deviations
@@ -546,13 +535,16 @@
     nIterations <- dim(scores)[1]
     nWaves <- dim(scores)[2]
     nParameters <- dim(scores)[3]
-    dfra <- rowSums(sapply(1:nWaves, function(i)
-                          sapply(1:nIterations, function(j, y, s)
-                                 outer(y[j, ], s[j, ]),
-                                 y=matrix(deviations[, i, ], ncol=nParameters),
-                                 s=matrix(scores[, i, ], ncol=nParameters))))
-    dim(dfra)<- c(nParameters, nParameters, nIterations)
-    dfra<- apply(dfra, c(1, 2), sum)
+    ## replaced nested applies by memory saving for loops after problems
+    dfra <- matrix(0, nrow=nParameters, ncol=nParameters)
+    for (i in 1:nIterations)
+    {
+        for (j in 1:nWaves)
+        {
+            dfra <- dfra + outer(scores[i, j, ], deviations[i, j, ])
+        }
+    }
+    dfra <- t(dfra)
     dfra<- dfra / nIterations
     tmp <- matrix(sapply(1 : nWaves, function(i)
                          outer(colMeans(deviations)[i,],
@@ -560,3 +552,21 @@
 
     dfra - matrix(rowSums(tmp), nrow=nParameters)
 }
+
+##@makeZsmall siena07 create a minimal version of z to pass between processors.
+makeZsmall <- function(z)
+{
+    zsmall <- NULL
+    zsmall$theta <- z$theta
+    zsmall$Deriv <- z$Deriv
+    zsmall$Phase <- z$Phase
+    zsmall$nit <- z$nit
+    zsmall$FinDiff.method <- z$FinDiff.method
+    zsmall$int2 <- z$int2
+    zsmall$cl <- z$cl
+    zsmall$maxlike <- z$maxlike
+    zsmall$cconditional <- z$cconditional
+    zsmall$condvar <- z$condvar
+    zsmall$pp <- z$pp
+    zsmall
+}

Modified: pkg/RSiena/R/phase2.r
===================================================================
--- pkg/RSiena/R/phase2.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/phase2.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -190,16 +190,8 @@
 {
     z$nit <- 0
     ac <- 0
-    zsmall <- NULL
-    zsmall$theta <- z$theta
-    zsmall$Deriv <- z$Deriv
-    zsmall$Phase <- z$Phase
-    zsmall$int2 <- z$int2
-    zsmall$FinDiff.method <- z$FinDiff.method
-    zsmall$cl <- z$cl
     xsmall <- NULL
-    zsmall$cconditional <- z$cconditional
-    zsmall$condvar <- z$condvar
+    zsmall <- makeZsmall(z)
     repeat
     {
         z$n <- z$n+1
@@ -266,7 +258,9 @@
         else
         {
             zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
-            fra <- rowMeans(sapply(zz, function(x) colSums(x$fra)- z$targets))
+            fra <- sapply(zz, function(x) colSums(x$fra)- z$targets)
+            dim(fra) <- c(z$pp, z$int)
+            fra <- rowMeans(fra)
             zz$OK <- sapply(zz, function(x) x$OK)
             if (!all(zz$OK))
             {
@@ -277,7 +271,7 @@
         if (x$maxlike)
         {
             z$phase2fras[subphase, ,z$nit] <- fra
-            z$rejectprops[subphase, , z$nit] <- zz$rejectprop
+         #   z$rejectprops[subphase, , z$nit] <- zz$rejectprop
         }
         if (z$nit %% 2 == 1)
         {
@@ -324,11 +318,11 @@
         z$theta <- zsmall$theta
         z$thav <- z$thav + zsmall$theta
         z$thavn <- z$thavn + 1
-        ## if (x$maxlike && !is.null(x$moreUpdates) && x$moreUpdates > 0)
-        ## {
-        ##    z <- doMoreUpdates(z, x, x$moreUpdates * subphase)
-        ##    zsmall$theta <- z$theta
-        ## }
+       # if (x$maxlike && !is.null(x$moreUpdates) && x$moreUpdates > 0)
+       # {
+       #     z <- doMoreUpdates(z, x, x$moreUpdates * subphase)
+       #     zsmall$theta <- z$theta
+       # }
         ##check for user interrupt
        ##   browser()
         CheckBreaks()

Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/phase3.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -64,17 +64,9 @@
         Report('Estimation of derivatives by the finite difference method.\n\n',outf)
     else
         Report('Estimation of derivatives by the LR method (type 1).\n\n', outf)
-    zsmall <- NULL
-    zsmall$FinDiff.method <- z$FinDiff.method
-    zsmall$theta <- z$theta
-    zsmall$Deriv <- z$Deriv
-    zsmall$Phase <- z$Phase
-    zsmall$nit <- z$nit
-    zsmall$int2 <- z$int2
-    zsmall$cl <- z$cl
+
     xsmall<- NULL
-    zsmall$cconditional <- z$cconditional
-    zsmall$condvar <- z$condvar
+    zsmall <- makeZsmall(z)
     if (!x$maxlike && !is.null(z$writefreq))
     {
         if (z$FinDiff.method)
@@ -195,6 +187,7 @@
                                  'phase-3 iterations.\n'), outf)
                     }
                     z$sf <- z$sf[1:nit,,drop=FALSE]
+                    z$sf2 <- z$sf2[1:nit,,,drop=FALSE]
                     z$ssc <- z$ssc[1:nit,,,drop=FALSE]
                     z$sdf <-z$sdf[1:nit,,,drop=FALSE]
                     break
@@ -209,6 +202,7 @@
         return(z)
     }
     z$Phase3nits <- nit
+    z$n3 <- nit
     z<- phase3.2(z,x)
     z
 }

Modified: pkg/RSiena/R/printDataReport.r
===================================================================
--- pkg/RSiena/R/printDataReport.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/printDataReport.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -131,7 +131,7 @@
             {
                 mdnet <- names(x$MaxDegree)[i]
                 Report(c("Dependent network variable", mdnet, ':\n'), outf)
-                if (attr(f, 'symmetric')[match(mdnet), attr(f, "netnames")])
+                if (attr(f, 'symmetric')[match(mdnet, attr(f, "netnames"))])
                 {
                     Report(c("All graphs are constrained to having degrees not",
                              "larger than", x$MaxDegree[i], '.\n'), outf)

Modified: pkg/RSiena/R/printInitialDescription.r
===================================================================
--- pkg/RSiena/R/printInitialDescription.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/printInitialDescription.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -99,27 +99,27 @@
             }
             averageOutDegree <- rep(NA, nData)
             for (group in 1:nData)
-                {
+            {
                 j <- match(netnames[net], names(data[[group]]$depvars))
                 if (is.na(j))
                     stop("network names not consistent")
                 depvar <- data[[group]]$depvars[[j]]
                 atts <- attributes(depvar)
                 averageOutDegree[group] <- atts$"averageOutDegree"
-                }
+            }
+            Report("\n", outf)
+            if (nData > 1)
+            {
+                Report("The average degrees are: ", outf)
+                Report(paste(names(data), round(averageOutDegree, 3),
+                             sep=': '), outf)
                 Report("\n", outf)
-            if (nData > 1)
-                {
-                    Report("The average degrees are: ", outf)
-                    Report(paste(names(data), round(averageOutDegree, 3),
-                                 sep=': '), outf)
-                    Report("\n", outf)
-                }
-                else
-                {
-                    Report(c("The average degree is",
-                             round(averageOutDegree, 3), "\n"), outf)
-                }
+            }
+            else
+            {
+                Report(c("The average degree is",
+                         round(averageOutDegree, 3), "\n"), outf)
+            }
             Report("\n\n", outf)
             Report(c(ifelse(gpatts$symmetric[net], "Edge", "Tie"),
                      "changes between subsequent observations:\n"), outf)

Modified: pkg/RSiena/R/robmon.r
===================================================================
--- pkg/RSiena/R/robmon.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/robmon.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -27,6 +27,7 @@
     z$maxrepeatsubphase <- 4
     z$gain <- x$firstg
     z$haveDfra <- FALSE
+    z$maxlike <- x$maxlike
     #######################################################
     ##do initial setup call of FRAN
     #######################################################
@@ -40,7 +41,7 @@
     #######################################################
     if (useCluster)
     {
-        if (!is.null(x$simstats0c) && !x$simstats0c)
+        if (!is.null(x$FRANname) && x$FRANname != "simstats0c")
         {
             stop("Multiple processors only for simstats0c at present")
         }

Modified: pkg/RSiena/R/siena01.r
===================================================================
--- pkg/RSiena/R/siena01.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/siena01.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -239,7 +239,7 @@
         ## browser()
         if (loadfilename == "")
         {
-            return()
+            return(FALSE)
         }
         modelName <<- basename(loadfilename)
         ipos <- max(c(0, gregexpr('.', modelName, fixed=TRUE)[[1]]))
@@ -249,26 +249,30 @@
         }
         session <<- sessionFromFile(loadfilename, tk=TRUE)
         procSession()
+        TRUE
     }
     ##@fromFileContFn internal siena01Gui
     fromFileContFn <- function()
     {
-        fromFileFn()
-        ## try to read in the project object
-        savedModelName <- paste(modelName, ".Rdata", sep='')
-        #browser()
-        if (inherits(resp <- try(load(savedModelName),
-                                 silent=TRUE), "try-error"))
+        OK <- fromFileFn()
+        if (OK)
         {
-            tkmessageBox(message="Unable to load saved model", icon="error")
+            ## try to read in the project object
+            savedModelName <- paste(modelName, ".Rdata", sep='')
+                                        #browser()
+            if (inherits(resp <- try(load(savedModelName),
+                                     silent=TRUE), "try-error"))
+            {
+                tkmessageBox(message="Unable to load saved model", icon="error")
+            }
+            else
+            {
+                mydata <<- mydata
+                mymodel <<- mymodel
+                myeff <<- myeff
+                sienaModelOptions()
+            }
         }
-        else
-        {
-            mydata <<- mydata
-            mymodel <<- mymodel
-            myeff <<- myeff
-            sienaModelOptions()
-        }
     }
     ##@helpFn internal siena01Gui
     helpFn <- function() ## display the manual
@@ -539,7 +543,6 @@
         {
             ##create mymodel
             mymodel <<- modelFromTcl()
-
             if (tclvalue(clustVar) == '0')
             {
                 nbrNodes <- 1
@@ -727,9 +730,54 @@
             if (savefilename != "")
                 save(estimAns, file=savefilename)
         }
-        ########################################
+        ##@screenFromModel internal siena01Gui update screen from saved model
+        screenFromModel <- function()
+        {
+            if (exists("mymodel"))
+            {
+                if (!is.na(mymodel$cconditional))
+                {
+                    if (mymodel$cconditional)
+                    {
+                        tclvalue(estimVar) <<-
+                            '1. conditional Method of Moments'
+                        if (ndepvars > 1)
+                        {
+                            tclvalue(condVar) <<- mymodel$condvarno
+                        }
+                    }
+                    else
+                    {
+                        tclvalue(estimVar) <<-
+                            '0. unconditional Method of Moments'
+                    }
+                }
+                tclvalue(gainVar) <<- mymodel$firstg
+                tclvalue(stdstartVar) <<- as.numeric(mymodel$useStdInits)
+                tclvalue(ph2spinVar) <<- mymodel$nsub
+                if (!is.null(mymodel$randomSeed) && mymodel$randomSeed != 0)
+                {
+                    tclvalue(rsVar) <<- '1'
+                    tclvalue(rsspinVar) <<- mymodel$randomSeed
+                    randomseedFn()
+                }
+                if (mymodel$FinDiff.method)
+                {
+                    tclvalue(derivVar) <<- '0. crude Monte Carlo'
+               }
+                tclvalue(ph3spinVar) <<- mymodel$n3
+                if (mymodel$MaxDegree != 0)
+                {
+                    for (i in 1:nMaxDegree)
+                    {
+                        maxdfVar[[i, 2]] <<- mymodel$MaxDegree[depvarnames[i]]
+                    }
+               }
+            }
+        }
+        ##*######################################
         ## start of sienaModelOptions function proper
-        ########################################
+        ##*#####################################
         resultsOpen <- FALSE
         if (inherits(mydata, "sienaGroup"))
         {
@@ -963,7 +1011,8 @@
         tkgrid(editbut, showbut,  resultsbut, saveresultsbut,
                 helpbut, padx=5, pady=5)
         tkgrid(effectsvarl, sticky="w")
-        ## make sure this window is top with a global grab, bu only for a second
+        screenFromModel()
+       ## make sure this window is top with a global grab, bu only for a second
         tcl('wm', 'attributes', tt, '-topmost', 1)
         Sys.sleep(0.1)
         tcl('wm', 'attributes', tt, '-topmost', 0)

Modified: pkg/RSiena/R/siena07.r
===================================================================
--- pkg/RSiena/R/siena07.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/siena07.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -1,7 +1,7 @@
 ##/*****************************************************************************
 ## * SIENA: Simulation Investigation for Empirical Network Analysis
 ## *
-## * Web: http://stat.gamma.rug.nl/siena.html
+## * Web: http://www.stats.ox.ac.uk/~snijders/siena
 ## *
 ## * File: siena07.r
 ## *

Modified: pkg/RSiena/R/sienaDataCreate.r
===================================================================
--- pkg/RSiena/R/sienaDataCreate.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/sienaDataCreate.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -473,7 +473,7 @@
         attr(depvars[[i]], 'distance') <- rep(FALSE, observations - 1)
         attr(depvars[[i]], 'vals') <- vector("list", observations)
         attr(depvars[[i]], 'nval') <- rep(NA, observations)
-        attr(depvars[[i]], 'noMissing') <- rep(NA, observations)
+        attr(depvars[[i]], 'noMissing') <- rep(0, observations)
         if (type == 'behavior')
         {
             attr(depvars[[i]], 'noMissing') <- FALSE
@@ -624,7 +624,7 @@
                     }
                     attr(depvars[[i]], "nval")[j] <-
                         sum(!is.na(mymat[row(mymat) != col(mymat)]))
-           }
+                }
                 ### need to exclude the structurals here
                 if (sparse)
                 {
@@ -656,12 +656,15 @@
                 degree <- (atts$netdims[1] - 1) * ones / atts$nval
                 missings <- 1 - atts$nval/ atts$netdims[1] /
                     (atts$netdims[1] - 1)
+                noMissing <- atts$netdims[1] * (atts$netdims[1] - 1) -
+                    atts$nval
                 attr(depvars[[i]], "ones") <- ones
                 attr(depvars[[i]], "density") <- density
                 attr(depvars[[i]], "degree") <- degree
                 attr(depvars[[i]], "averageOutDegree") <- mean(degree)
                 attr(depvars[[i]], "averageInDegree") <- mean(degree)
                 attr(depvars[[i]], "missings") <- missings
+                attr(depvars[[i]], "noMissing") <- noMissing
             }
             else #type=='bipartite' not sure what we need here,
                 ## but include diagonal
@@ -734,12 +737,14 @@
                 density <- ones / atts$nval
                 degree <- (atts$netdims[2]) * ones / atts$nval
                 missings <- 1 - atts$nval/ atts$netdims[1] /
-                    (atts$netdims[2])
+                    atts$netdims[2]
+                noMissing <- atts$netdims[1] * atts$netdims[2] - atts$nval
                 attr(depvars[[i]], "ones") <- ones
                 attr(depvars[[i]], "density") <- density
                 attr(depvars[[i]], "degree") <- degree
                 attr(depvars[[i]], "averageOutDegree") <- mean(degree)
                 attr(depvars[[i]], "missings") <- missings
+                attr(depvars[[i]], "noMissing") <- noMissing
            }
         }
         attr(depvars[[i]], 'name') <- names(depvars)[i]
@@ -754,11 +759,7 @@
     z$dycCovars <- dycCovars
     z$dyvCovars <- dyvCovars
     z$compositionChange <- compositionChange
- #   types <- sapply(z$depvars, function(x)attr(x, "type"))
- #   if (sum(types != "behavior" ) > 1)
- #   {
-        z <- checkConstraints(z)
- #   }
+    z <- checkConstraints(z)
     class(z) <- 'siena'
     z
 }
@@ -767,36 +768,42 @@
 {
     types <- sapply(z$depvars, function(x)attr(x, "type"))
     sparse <- sapply(z$depvars, function(x)attr(x, "sparse"))
-    nodeSets <- sapply(z$depvars, function(x)attr(x, "nodeSet"))
+    nodeSets <- lapply(z$depvars, function(x)attr(x, "nodeSet"))
     nNets <- length(z$depvars)
+
     pairsOfNets <- as.matrix(expand.grid(names(z$depvars), names(z$depvars)))
-    ## maybe remove some as don't want pairs with self, but may want all there
-    ##pairsOfNets <- pairsOfNets[pairsOfNets[, 1] != pairsOfNets[, 2], ]
-    pairsNames <- paste(pairsOfNets[, 1], pairsOfNets[,2], sep=",")
+    pairsNames <- paste(pairsOfNets[, 1], pairsOfNets[, 2], sep=",")
 
     higher <- namedVector(FALSE, pairsNames )
     atLeastOne <- namedVector(FALSE, pairsNames )
     disjoint <- namedVector(FALSE, pairsNames )
 
-    ## identify any nets which may relate
+    ## identify any nets which may relate. These are those that
+    ## share a node set and type.
     relates <- data.frame(name=names(z$depvars), type=types,
                           nodeSets=sapply(nodeSets, paste, collapse=","),
                           tn=paste(types, sapply(nodeSets, paste,
                           collapse=",")) , stringsAsFactors=FALSE)
+
+    ## just check we have only one row per network
+    if (nrow(relates) != nNets)
+    {
+        stop("Error in checkConstraints")
+    }
+    ## find the ones that do occur more than once
     use <- relates$tn %in% relates$tn[duplicated(relates$tn)]
+
+    ## create a working list to be compared, with names for ease of access
     nets <- namedVector(NA, names(z$depvars), listType=TRUE)
+
     for (net in names(z$depvars)[use])
     {
         if (types[[net]] != "behavior")
         {
             nets[[net]] <- z$depvars[[net]]
-        ##    nets[[net]] <- replaceMissingsAndStructurals(z$depvars[[net]])
         }
     }
 
-   ## relSplits <- split(relates, relates$tn)
-   ## relSplits <- relSplits[sapply(relSplits, nrow) > 1]
-
     for (i in 1:nrow(pairsOfNets))
     {
         if (pairsOfNets[i, 1] != pairsOfNets[i, 2])
@@ -809,7 +816,7 @@
             nodes1 <- relates[net1, "nodeSets"]
             nodes2 <- relates[net2, "nodeSets"]
 
-            if (type1 == type2 && type1 != "behavior" & nodes1 == nodes2)
+            if (type1 == type2 && type1 != "behavior" && nodes1 == nodes2)
             {
                 higher[i] <- TRUE
                 disjoint[i] <- TRUE
@@ -856,7 +863,6 @@
             }
         }
     }
-
     attr(z, "higher") <- higher
     attr(z, "disjoint") <- disjoint
     attr(z, "atLeastOne") <- atLeastOne
@@ -875,14 +881,6 @@
         rvals <- range(vals,na.rm=TRUE)
     }
     rvals1 <- rvals[2] - rvals[1]
-    ## use of outer creates a potentially huge matrix unnecessarily
-    ## tmp3 <- apply(vals,2,function(x)
-    ##          {
-    ##              y <- 1 - abs(outer(x, x, "-")) / rvals1
-    ##              diag(y) <- NA
-    ##              y
-    ##          })
-    ## tmp3 <- array(tmp3, dim=c(n, n, observations))
     tmp <- apply(vals, 2, function(v){
     sapply(1: length(v), function(x, y, r){
         z <- y
@@ -1157,7 +1155,7 @@
     attr(group, "dyvCovarMean") <- dyvCovarMean
     group
 }
-##@namedVector DataCreate
+##@namedVector DataCreate Utility to create a vector with appropriate names
 namedVector <- function(vectorValue, vectorNames, listType=FALSE)
 {
     if (listType)
@@ -1255,6 +1253,7 @@
     cvnodeSets <- namedVector(NA, vCovars)
     dycnodeSets <- namedVector(NA, dycCovars, listType=TRUE)
     dyvnodeSets <- namedVector(NA, dyvCovars, listType=TRUE)
+    totalMissings <- namedVector(0, netnames)
     observations <- 0
     periodNos <- rep(NA, 2)
     groupPeriods <- namedVector(NA, names(objlist))
@@ -1314,6 +1313,8 @@
             {
                 structural[netnamesub] <- TRUE
             }
+            totalMissings[netnamesub] <- totalMissings[netnamesub] +
+                sum(attribs[["noMissing"]])
         }
         thisHigher <- attr(objlist[[i]], "higher")
         thisDisjoint <- attr(objlist[[i]], "disjoint")
@@ -1472,6 +1473,7 @@
     attr(group, 'netnames') <- netnames
     attr(group, 'symmetric') <- symmetric
     attr(group, 'structural') <- structural
+    attr(group, "totalMissings") <- totalMissings
     attr(group, 'allUpOnly') <- allUpOnly
     attr(group, 'allDownOnly') <- allDownOnly
     attr(group, 'anyUpOnly') <- anyUpOnly

Modified: pkg/RSiena/R/sienaDataCreateFromSession.r
===================================================================
--- pkg/RSiena/R/sienaDataCreateFromSession.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/sienaDataCreateFromSession.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -525,7 +525,7 @@
                      ##  namefiles[[1]][grep(miss, namefiles[[1]])] <-  NA
                        miss <- namesession$MissingValues
                         miss <- strsplit(miss, " ")[[1]]
-                       namefiles[[1]][namefiles[[1]] %in% miss] <-  NA
+                      namefiles[[1]][namefiles[[1]] %in% miss] <-  NA
                        if (namesession[1, "ActorSet"] == "Actors")
                        {
                            namesession[1, "ActorSet"]<- "Actors Actors"

Modified: pkg/RSiena/R/sienaModelCreate.r
===================================================================
--- pkg/RSiena/R/sienaModelCreate.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/sienaModelCreate.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -11,12 +11,13 @@
 
 ##@sienaModelCreate DataCreate
 sienaModelCreate <-
-    function(fn=simstats0c,
-             usesimstats0c=deparse(substitute(fn)) == "simstats0c",
+    function(fn,
              projname="Siena", MaxDegree=0, useStdInits=FALSE,
              n3=1000, nsub=4, maxlike=FALSE, diag=!maxlike,
              condvarno=0, condname='',
-             firstg=0.2, cond=NA, findiff=FALSE,  seed=NULL)
+             firstg=0.2, cond=NA, findiff=FALSE,  seed=NULL,
+             pridg=0.05, prcdg=0.05, prper=0.3, pripr=0.25, prdpr=0.25,
+             prrms=0.1)
 {
     model <- NULL
     model$projname <- projname
@@ -26,9 +27,36 @@
     model$firstg <- firstg
     model$maxrat <- 1.0
     model$maxmaxrat <- 10.0
-   # model$FRAN <- fn
+    model$maxlike <-  maxlike
     model$FRANname <- deparse(substitute(fn))
-    model$maxlike <-  maxlike
+    if (maxlike)
+    {
+        if (missing(fn))
+        {
+            model$FRANname <- "maxlikec"
+        }
+        if (is.na(cond))
+        {
+            cond <- FALSE
+        }
+        if (cond)
+        {
+            stop("Conditional estimation is not possible with",
+                  "maximum likelihood estimation")
+        }
+        if (findiff)
+        {
+            stop("Finite differences estimation of derivatives",
+                 "is not possible with maximum likelihood estimation")
+        }
+    }
+    else
+    {
+        if (missing(fn))
+        {
+            model$FRANname <- "simstats0c"
+        }
+    }
     model$cconditional <- cond
     model$condvarno <-  condvarno
     model$condname <- condname
@@ -38,7 +66,12 @@
     model$ModelType <- 1
     model$MaxDegree <- MaxDegree
     model$randomSeed <- seed
-    model$simstats0c <- usesimstats0c
+    model$pridg <- pridg
+    model$prcdg <- prcdg
+    model$prper <- prper
+    model$pripr <- pripr
+    model$prdpr <- prdpr
+    model$prrms <- prrms
     class(model) <- "sienaModel"
     model
 }

Modified: pkg/RSiena/R/sienaprint.r
===================================================================
--- pkg/RSiena/R/sienaprint.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/sienaprint.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -253,7 +253,7 @@
             mydf[1:nn, 'row'] <- nnstr
             mydf[1:nn, 'value'] <- x$rate
             mydf[1:nn, 'se'] <- x$vrate
-            if (x$f$nDepvars == 1 || is.null(x$f$nDepvars))
+            if (is.null(x$f[[1]]$nDepvars) || x$f[[1]]$nDepvars == 1)
             {
                 mydf[1:nn, 'text'] <- paste('Rate parameter period', 1:nn)
             }

Modified: pkg/RSiena/R/simstatsc.r
===================================================================
--- pkg/RSiena/R/simstatsc.r	2010-04-12 16:01:03 UTC (rev 76)
+++ pkg/RSiena/R/simstatsc.r	2010-04-12 16:38:28 UTC (rev 77)
@@ -6,7 +6,7 @@
 # * File: simstatsc.r
 # *
 # * Description: This module contains the code for simulating the process,
-# * communicating with C++.
+# * communicating with C++. Only subsidiary routines used for maximum likelihood
 # *****************************************************************************/
 ##@simstats0c siena07 Simulation Module
 simstats0c <-function(z, x, INIT=FALSE, TERM=FALSE, initC=FALSE, data=NULL,
@@ -52,15 +52,12 @@
         dimnames(z$dfra)[[1]] <- as.list(z$requestedEffects$shortName)
         return(z)
     }
-    ## iteration entry points
+    ######################################################################
+    ## iteration entry point: not for ML
+    ######################################################################
+    ## retrieve stored information
     f <- FRANstore()
     ## browser()
-    ## cat(f$randomseed2, f$storedseed, '\n')
-    ## if (fromFiniteDiff)
-    ## {
-    ##  cat('before\n')
-    ##  print(f$seeds)
-    ## }
     if (fromFiniteDiff || z$Phase == 2)
     {
         returnDeps <- FALSE
@@ -94,17 +91,21 @@
         }
        ## cat(randomseed2, '\n')
     }
+    ## create a grid of periods with group names in case want to parallelize
+    ## using this
     groupPeriods <- attr(f, "groupPeriods")
     callGrid <- cbind(rep(1:f$nGroup, groupPeriods - 1),
                       as.vector(unlist(sapply(groupPeriods - 1,
                                               function(x) 1:x))))
+    ## z$int2 is the number of processors if iterating by period, so 1 means
+    ## we are not
     if (z$int2==1 || nrow(callGrid) == 1)
     {
       #   cat("theta", z$theta, "\n")
-       ans <- .Call('model', PACKAGE=pkgname,
-                     z$Deriv, f$pData, seeds,
-                     fromFiniteDiff, f$pModel, f$myeffects, z$theta,
-                     randomseed2, returnDeps, z$FinDiff.method, !is.null(z$cl))
+            ans <- .Call('model', PACKAGE=pkgname, z$Deriv, f$pData, seeds,
+                         fromFiniteDiff, f$pModel, f$myeffects, z$theta,
+                         randomseed2, returnDeps, z$FinDiff.method,
+                         !is.null(z$cl))
     }
     else
     {
@@ -125,19 +126,21 @@
         ans[[5]] <- NULL # randomseed not sensible here
         fff <- lapply(anss, "[[", 6)
         fff <- split(fff, callGrid[, 1])
-        ans[[6]] <- lapply(fff, function(x){
-            lapply(1:length(f$depNames), function(x, z) lapply(z, "[[", x), z=x)
-        })
+        ans[[6]] <-
+            lapply(fff, function(x)
+               {
+                   lapply(1:length(f$depNames), function(x, z)
+                          lapply(z, "[[", x), z=x)
+               }
+                   )
     }
     ## browser()
     if (!fromFiniteDiff)
     {
         if (z$FinDiff.method)
             f$seeds <- ans[[3]]
-        ## cat('after\n')
-        ## print(f$seeds)
     }
-    if (z$Deriv)
+    if (z$Deriv )
     {
         sc <- t(ans[[2]])
     }
@@ -169,24 +172,30 @@
             periodNo <- periodNos[length(periodNos)] + 2
        }
     }
-    ##   cat('fra', fra, '\n')
-    ##    cat(f$randomseed2, f$storedseed, '\n')
+   # browser()
    list(sc = sc, fra = fra, ntim0 = ntim, feasible = TRUE, OK = TRUE,
          sims=sims, f$seeds)
 }
 doModel <- function(x, Deriv, seeds, fromFiniteDiff, theta, randomseed2,
-                    returnDeps, FinDiff.method, useStreams)
+                    returnDeps, FinDiff.method, useStreams, maxlike)
 {
-    ## cat('gothere\n')
-    ##     cat("theta", theta, "\n")
     f <- FRANstore()
-    ## cat('stream', .lec.GetStreams(),'\n')
-    ## print(.lec.GetState(.lec.GetStreams()))
     seeds <- seeds[[x[1]]][[x[2]]]
-    .Call("modelPeriod", PACKAGE=pkgname, Deriv, f$pData, seeds,
-                         fromFiniteDiff, f$pModel, f$myeffects, theta,
-                         randomseed2, returnDeps, FinDiff.method, useStreams,
-          as.integer(x[1]), as.integer(x[2]))
+    if (!maxlike)
+    {
+        .Call("modelPeriod", PACKAGE=pkgname, Deriv, f$pData, seeds,
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/rsiena -r 77


More information about the Rsiena-commits mailing list