[Rsiena-commits] r180 - in pkg: RSiena RSiena/R RSiena/man RSiena/src/model RSiena/src/model/effects RSienaTest RSienaTest/R RSienaTest/man RSienaTest/src/model

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Nov 4 13:33:30 CET 2011


Author: ripleyrm
Date: 2011-11-04 13:33:30 +0100 (Fri, 04 Nov 2011)
New Revision: 180

Modified:
   pkg/RSiena/DESCRIPTION
   pkg/RSiena/R/phase1.r
   pkg/RSiena/R/phase3.r
   pkg/RSiena/changeLog
   pkg/RSiena/man/RSiena-package.Rd
   pkg/RSiena/src/model/StatisticCalculator.cpp
   pkg/RSiena/src/model/effects/BehaviorInteractionEffect.cpp
   pkg/RSienaTest/DESCRIPTION
   pkg/RSienaTest/R/phase1.r
   pkg/RSienaTest/R/phase3.r
   pkg/RSienaTest/changeLog
   pkg/RSienaTest/man/RSiena-package.Rd
   pkg/RSienaTest/src/model/StatisticCalculator.cpp
Log:
Created a function to do block of iterations for phases 1 or 3. Removed functions to do just one and to update the stores as slow. Added if statements to StatisticCalculator as a quick speed up. 

Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION	2011-10-28 12:43:23 UTC (rev 179)
+++ pkg/RSiena/DESCRIPTION	2011-11-04 12:33:30 UTC (rev 180)
@@ -1,8 +1,8 @@
 Package: RSiena
 Type: Package
 Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.0.12.179
-Date: 2011-10-28
+Version: 1.0.12.180
+Date: 2011-11-04
 Author: Various
 Depends: R (>= 2.10.0)
 Imports: Matrix

Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r	2011-10-28 12:43:23 UTC (rev 179)
+++ pkg/RSiena/R/phase1.r	2011-11-04 12:33:30 UTC (rev 180)
@@ -5,16 +5,16 @@
 ## *
 ## * File: phase1.r
 ## *
-## * Description: This module contains the functions phase1.1, doPhase1it,
+## * Description: Phase 1 processing is controlled by the functions phase1.1
 ## * and phase1.2. phase1.1 does the first 10 iterations and checks that
 ## * finite differencing is going OK. If not, alters epsilon and forces a
 ## * restart. Phase 1.2 does the rest of the iterations and then calculates
-## * the derivative estimate. doPhase1it does one iteration, is called by
-## * phase1.1 and phase1.2.
+## * the derivative estimate. Both call doPhase1or3iterations which does a block
+## * of iterations, to save copying large objects around.
 ## ****************************************************************************/
 ##args: x model object (readonly), z control object
 ##
-##@phase1.1 siena07 Do first 10 iterations (before check if using finite differences)
+##@phase1.1 siena07 Do first 10 iters (before check if using finite differences)
 phase1.1 <- function(z, x, ...)
 {
     ## initialise phase 1
@@ -55,45 +55,11 @@
         nits6 <- min(nits[nits >= 6 ])
     }
     else
+	{
         nits6 <- -1
-    for (nit in nits)
-    {
-         z$nit <- nit
-         if (length(nits) > 1 && nit == nits[2])
-         {
-             time1 <- proc.time()['elapsed']
-             if (x$checktime)
-             {
-                 z$ctime <- time1
-             }
-         }
-         if (nit == nits6)
-         {
-             time1 <- proc.time()['elapsed']-time1
-             if (time1 > 1e-5)
-             {
-                 z$writefreq <- round(20.0/time1)
-             }
-             else
-             {
-                 z$writefreq <- 5
-             }
-             if (is.batch())
-             {
-                 z$writefreq <- 10 * z$writefreq
-             }
-             z$writefreq <- roundfreq(z$writefreq)
-         }
-         ## #######################################
-         ## # do iteration
-         ## #######################################
-         z <- doPhase1it(z, x, zsmall=zsmall, xsmall=xsmall, ...)
-         ## #######################################
-         ## #
-         ## #######################################
-         if (!z$OK || UserRestartFlag() || UserInterruptFlag())
-             return(z)
-     }
+	}
+	### call the subroutine to do the iterations
+	z <- doPhase1or3Iterations(1, z, x, zsmall, xsmall, nits, nits6)
 
     if (z$FinDiff.method)
     {
@@ -153,66 +119,7 @@
     z
 }
 
-##@doPhase1it siena07 does 1 iteration in Phase 1
-doPhase1it<- function(z, x, zsmall, xsmall, ...)
-{
-    DisplayIteration(z)
-    if (z$int == 1)
-    {
-        zz <- x$FRAN(zsmall, xsmall)
-        if (!zz$OK)
-        {
-            z$OK <- zz$OK
-            z$zz <- zz
-            return(z)
-        }
-        z$n <- z$n + 1
-        z$phase1Its <- z$phase1Its + 1
-    }
-    else
-    {
-        zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
-        z$n <- z$n + z$int
-        z$phase1Its <- z$phase1Its + z$int
-      #  browser()
-    }
-	z <- updateSiena07stores(z, zz, x)
-    CheckBreaks()
-    if (UserInterruptFlag() || UserRestartFlag())
-    {
-        return(z)
-    }
-    val <- getProgressBar(z$pb)
-    if (z$FinDiff.method)
-    {
-        val <- val + z$pp + 1
 
-    }
-    else
-    {
-        val <- val + 1
-    }
-    z$pb <- setProgressBar(z$pb, val)
-    progress <- val / z$pb$pbmax * 100
-    if (z$nit <= 5 || z$nit %% z$writefreq == 0 || z$nit %%5 == 0 ||
-        x$maxlike || z$FinDiff.method ||
-        (z$int > 1 && z$nit %% z$writefreq < z$int))
-    {
-      #  Report(c('Phase', z$Phase, 'Iteration ', z$nit, '\n'))
-        if (is.batch())
-		{
-            Report(c('Phase ', z$Phase, ' Iteration ', z$nit, ' Progress: ',
-                     round(progress), '%\n'), sep='')
-		}
-        else
-        {
-            DisplayTheta(z)
-            DisplayDeviations(z, z$sf[z$nit, ])
-        }
-    }
-    #browser()
-    z
-}
 ##@phase1.2 siena07 Do rest of phase 1 iterations
 phase1.2 <- function(z, x, ...)
 {
@@ -223,23 +130,7 @@
     if (z$n1 > z$phase1Its)
     {
         nits <- seq((z$phase1Its+1), z$n1, int)
-        for (nit in nits)
-        {
-            if (is.null(z$ctime))
-            {
-                time1 <- proc.time()['elapsed']
-                if (x$checktime)
-                {
-                    z$ctime <- time1
-                }
-            }
-            z$nit <- nit
-            z <- doPhase1it(z, x, zsmall=zsmall, xsmall=xsmall, ...)
-            if (!z$OK || UserInterruptFlag() || UserRestartFlag())
-            {
-                return(z)
-            }
-        }
+		z <- doPhase1or3Iterations(1, z, x, zsmall, xsmall, nits)
     }
     z$timePhase1 <- (proc.time()['elapsed'] - z$ctime) / (z$nit - 1)
     if (x$checktime  && !is.na(z$timePhase1))
@@ -252,7 +143,7 @@
     Report('after phase 1:\n', cf)
     PrtOutMat(format(as.matrix(z$mnfra), width = 15, nsmall = 6), cf)
     z <- CalculateDerivative(z, x)
-  ##browser()
+	##browser()
     if (!z$OK || z$DerivativeProblem) ##longer phase 1 or use finite differences
     {
         return(z)
@@ -567,88 +458,3 @@
     z$sims <- vector("list", nIterations)
 	z
 }
-##@updateSiena07Stores siena07 store data in phase 1 and 3
-updateSiena07stores <- function(z, zz, x)
-{
-	int <- z$int
-	if (int == 1)
-    {
-        fra <- colSums(zz$fra)
-        fra <- fra - z$targets
-        fra2 <- zz$fra
-        z$sf[z$nit, ] <- fra
-        z$sf2[z$nit, , ] <- zz$fra
-        z$sims[[z$nit]] <- zz$sims
-        z$chain[[z$nit]] <- zz$chain
-        fra <- fra + z$targets
-    }
-    else
-    {
-        for (i in 1:int)
-        {
-            fra <- colSums(zz[[i]]$fra)
-            fra <- fra - z$targets
-            z$sf[z$nit + (i - 1), ] <- fra
-            z$sf2[z$nit + (i - 1), , ] <- zz[[i]]$fra
-            z$sims[[z$nit + (i - 1)]] <- zz[[i]]$sims
-        }
-        fra2 <- t(sapply(zz, function(x)x$fra))
-        dim(fra2) <- c(int, nrow(zz[[1]]$fra), z$pp)
-        fra <- t(sapply(zz, function(x) colSums(x$fra)))
-    }
-    if (x$maxlike)
-    {
-        z$sdf[[z$nit]] <- zz$dff
-        z$sdf2[[z$nit]] <- zz$dff2
-        z$accepts[z$nit, , ] <- zz$accepts
-        z$rejects[z$nit, , ] <- zz$rejects
-		z$aborts[z$nit, , ] <- zz$aborts
-	}
-    else if (z$FinDiff.method)
-    {
-        z <- FiniteDifferences(z, x, fra, fra2)
-		for (i in 0:(z$int - 1))
-		{
-			z$sdf[[z$nit + i]] <- z$sdf0[i + 1, , ]
-			if (z$byWave)
-			{
-				z$sdf2[[z$nit + i]] <- z$sdf02[i + 1, , , ]
-			}
-		}
-    }
-    else
-    {
-        if (int==1)
-        {
-            if (!is.null(zz[['sc']]))
-			{
-                z$ssc[z$nit , ,] <- zz$sc
-			}
-        }
-        else
-        {
-                for (i in 1:int)
-                {
-                    if (!is.null(zz[[i]][['sc']]))
-					{
-                        z$ssc[z$nit + (i - 1), , ] <- zz[[i]]$sc
-					}
-                }
-            }
-    }
-    if ((!x$maxlike) && z$cconditional && z$Phase == 3)
-    {
-        if (int == 1)
-        {
-            z$ntim[z$nit, ] <- zz$ntim0
-        }
-        else
-        {
-            for (i in 1:int)
-            {
-				z$ntim[z$nit + (i-1), ] <- zz[[i]]$ntim0
-            }
-        }
-    }
-	z
-}

Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r	2011-10-28 12:43:23 UTC (rev 179)
+++ pkg/RSiena/R/phase3.r	2011-11-04 12:33:30 UTC (rev 180)
@@ -33,6 +33,7 @@
         endNit <- x$n3
     }
 	z$n3 <- endNit
+    z$Phase3nits <- endNit
 
     ## revert to original requested method for phase 3 unless symmetric
     if (z$FinDiff.method && !x$FinDiff.method &&
@@ -72,155 +73,12 @@
     {
         z$writefreq <- 10
     }
-    for (nit in nits)
-    {
-        z$nit <- nit
-        if (is.null(writefreq))
-        {
-            if (nit == nits[2])
-            {
-                z$time1 <- proc.time()[[3]]
-            }
-            else if (nit == nits11)
-            {
-                time1 <- proc.time()[[3]] - z$time1
-                if (time1 > 1e-5)
-                {
-                    z$writefreq <- max(1, round(20.0 / time1))
-                }
-                else
-                {
-                    z$writefreq <- 20
-                }
-                if (is.batch())
-                {
-                    z$writefreq <-  z$writefreq * 2 ##compensation for it
-                    ## running faster with no tcl/tk
-                }
-                z$writefreq <- roundfreq(z$writefreq)
-                writefreq <- z$writefreq
-            }
-        }
-        if (nit <= 5 || nit == 10 || (int==1 && nit %% z$writefreq == 0 ) ||
-            (int > 1 && (z$writefreq + 1) < z$n3%/%int &&
-             nit %in% nits[seq(z$writefreq + 1, x$n3 %/% int,
-                                          z$writefreq)]))
-        {
-            if (!is.batch())
-            {
-                DisplayIteration(z)
-              #  DisplayTheta(z)
-                DisplayDeviations(z, z$sf[nit-1,])
-            }
-           ## else
-           ## {
-           ##     Report(c('Phase ', z$Phase,' Iteration ',nit,'\n'))
-           ## }
-            if (nit %% z$writefreq == 0 || (int > 1 &&
-                       nit %% z$writefreq == 1) )
-            {
-                increment <- ifelse(nit <= 5, int,
-                                    ifelse(nit <= 10, 5, z$writefreq * int))
-                val<- getProgressBar(z$pb)
-                if (z$FinDiff.method)
-                    val <- val + increment * (z$pp + 1)
-                else
-                    val <- val + increment
-                ## sort out progress bar
-                z$pb <- setProgressBar(z$pb, val)
-                if (is.batch())
-                    Report(c("Phase ", z$Phase, " Iteration ", nit,
-                             " Progress ",
-                             round(val / z$pb$pbmax * 100), "%\n"), sep='')
-            }
-        }
-###############################
-        z <- doPhase3it(z, x, nit, zsmall=zsmall, xsmall=xsmall, ...)
-##############################
-        if (!is.batch())
-        {
-            if (nit < 10)
-                Report(c("  ", nit, " ", format(z$sf[nit,], width=10,
-                                             digits=4, nsmall=4), "\n"), cf)
-            if (nit >= 10)
-            {
-                CheckBreaks()
-                ##    if (nit==10) set up stopkey hint
-				##   Early termination of estimation
-                if( UserInterruptFlag())
-                {
-                   Report(c("The user asked for an early stop of the ",
-                             "algorithm during phase 3.\n"), outf)
-                    z$Phase3Interrupt <- TRUE
-                    if (nit < 500)
-                    {
-                        if (EarlyEndPhase2Flag())
-						{
-                            Report('This implies that ', outf)
-						}
-                        else
-						{
-                            Report(c("This implies that the estimates are as",
-								   "usual,\nbut the "), outf)
-						}
-                        Report(c('diagnostic checks, covariance matrices and ',
-                                 't-values \nare less reliable, because they ',
-                                 'are now based on only', nit,
-                                 'phase-3 iterations.\n\n'), outf)
-                    }
-                    z$sf <- z$sf[1:nit, , drop=FALSE]
-                    z$sf2 <- z$sf2[1:nit, , , drop=FALSE]
-                    if (!x$maxlike)
-                    {
-                        z$ssc <- z$ssc[1:nit, , , drop=FALSE]
-                    }
-                    else
-                    {
-                        z$sdf <-z$sdf[1:nit, , , drop=FALSE]
-                        z$sdf2 <-z$sdf2[1:nit, , , ,drop=FALSE]
-                    }
-                    z$sims <-z$sims[1:nit]
-                    endNit <- nit
-                    break
-                }
-                if (UserRestartFlag())
-                    break
-            }
-        }
-    }
-    if (!z$OK || UserRestartFlag())
-    {
-        return(z)
-    }
-    z$Phase3nits <- endNit
-    z$n3 <- endNit
+	z <- doPhase1or3Iterations(3, z, x, zsmall, xsmall, nits, 0, nits11,
+							   writefreq)
     z <- phase3.2(z,x)
     z
 }
 
-##@doPhase3it siena07 Does one iteration in phase 3
-doPhase3it<- function(z, x, nit, zsmall, xsmall, ...)
-{
-    if (z$int == 1)
-    {
-        zz <- x$FRAN(zsmall, xsmall)
-        if (!zz$OK)
-        {
-            z$OK <- zz$OK
-            z$zz <- zz
-            return(z)
-        }
-        z$n <- z$n + 1
-    }
-    else
-    {
-		zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
-        z$n <- z$n + z$int
- 	}
-	z <- updateSiena07stores(z, zz, x)
-    z
-}
-
 ##@phase3.2 siena07 Processing at end of phase 3
 phase3.2 <- function(z, x, ...)
 {
@@ -498,3 +356,320 @@
     }
     z
 }
+
+doPhase1or3Iterations <- function(phase, z, x, zsmall, xsmall, nits, nits6=0,
+								  nits11=0, writefreq)
+{
+	int <- z$int
+	for (nit in nits)
+	{
+		z$nit <- nit
+		## initial steps differ between phases
+		if (phase == 1)
+		{
+			if (length(nits) > 1 && nit == nits[2])
+			{
+				time1 <- proc.time()['elapsed']
+				if (x$checktime)
+				{
+					z$ctime <- time1
+				}
+			}
+			if (nit == nits6)
+			{
+				time1 <- proc.time()['elapsed']-time1
+				if (time1 > 1e-5)
+				{
+					z$writefreq <- round(20.0/time1)
+				}
+				else
+				{
+					z$writefreq <- 5
+				}
+				if (is.batch())
+				{
+					z$writefreq <- 10 * z$writefreq
+				}
+				z$writefreq <- roundfreq(z$writefreq)
+			}
+			if (nit > nits6 && is.null(z$ctime))
+			{
+				time1 <- proc.time()['elapsed']
+				if (x$checktime)
+				{
+					z$ctime <- time1
+				}
+			}
+			DisplayIteration(z)
+		}
+		else ## phase 3
+		{
+			if (is.null(writefreq))
+			{
+				if (nit == nits[2])
+				{
+					z$time1 <- proc.time()[[3]]
+				}
+				else if (nit == nits11)
+				{
+					time1 <- proc.time()[[3]] - z$time1
+					if (time1 > 1e-5)
+					{
+						z$writefreq <- max(1, round(20.0 / time1))
+					}
+					else
+					{
+						z$writefreq <- 20
+					}
+					if (is.batch())
+					{
+						z$writefreq <-  z$writefreq * 2 ##compensation for it
+						## running faster with no tcl/tk
+					}
+					z$writefreq <- roundfreq(z$writefreq)
+					writefreq <- z$writefreq
+				}
+			}
+			if (nit <= 5 || nit == 10 || (int==1 && nit %% z$writefreq == 0 ) ||
+				(int > 1 && (z$writefreq + 1) < z$n3%/%int &&
+				 nit %in% nits[seq(z$writefreq + 1, x$n3 %/% int,
+								   z$writefreq)]))
+			{
+				if (!is.batch())
+				{
+					DisplayIteration(z)
+					##  DisplayTheta(z)
+					DisplayDeviations(z, z$sf[nit-1,])
+				}
+				## else
+				## {
+				##     Report(c('Phase ', z$Phase,' Iteration ',nit,'\n'))
+				## }
+				if (nit %% z$writefreq == 0 || (int > 1 &&
+						   nit %% z$writefreq == 1))
+				{
+					increment <- ifelse(nit <= 5, int,
+										ifelse(nit <= 10, 5, z$writefreq * int))
+					val<- getProgressBar(z$pb)
+					if (z$FinDiff.method)
+						val <- val + increment * (z$pp + 1)
+					else
+						val <- val + increment
+					## sort out progress bar
+					z$pb <- setProgressBar(z$pb, val)
+					if (is.batch())
+						Report(c("Phase ", z$Phase, " Iteration ", nit,
+								 " Progress ",
+								 round(val / z$pb$pbmax * 100), "%\n"), sep='')
+				}
+			}
+		}
+		## iteration proper
+		if (z$int == 1)
+		{
+			zz <- x$FRAN(zsmall, xsmall)
+			if (!zz$OK)
+			{
+				z$OK <- zz$OK
+				z$zz <- zz
+				return(z)
+			}
+			z$n <- z$n + 1
+			if (phase == 1)
+			{
+				z$phase1Its <- z$phase1Its + 1
+			}
+		}
+		else
+		{
+			zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
+			z$n <- z$n + z$int
+			if (phase == 1)
+			{
+				z$phase1Its <- z$phase1Its + z$int
+			}
+		}
+		## now update the stores
+		if (z$int == 1)
+		{
+			fra <- colSums(zz$fra)
+			fra <- fra - z$targets
+			fra2 <- zz$fra
+			z$sf[z$nit, ] <- fra
+			z$sf2[z$nit, , ] <- zz$fra
+			z$sims[[z$nit]] <- zz$sims
+			z$chain[[z$nit]] <- zz$chain
+			fra <- fra + z$targets
+		}
+		else
+		{
+			for (i in 1:int)
+			{
+				fra <- colSums(zz[[i]]$fra)
+				fra <- fra - z$targets
+				z$sf[z$nit + (i - 1), ] <- fra
+				z$sf2[z$nit + (i - 1), , ] <- zz[[i]]$fra
+				z$sims[[z$nit + (i - 1)]] <- zz[[i]]$sims
+			}
+			fra2 <- t(sapply(zz, function(x)x$fra))
+			dim(fra2) <- c(int, nrow(zz[[1]]$fra), z$pp)
+			fra <- t(sapply(zz, function(x) colSums(x$fra)))
+		}
+		if (x$maxlike)
+		{
+			z$sdf[[z$nit]] <- zz$dff
+			z$sdf2[[z$nit]] <- zz$dff2
+			z$accepts[z$nit, , ] <- zz$accepts
+			z$rejects[z$nit, , ] <- zz$rejects
+			z$aborts[z$nit, , ] <- zz$aborts
+		}
+		else if (z$FinDiff.method)
+		{
+			z <- FiniteDifferences(z, x, fra, fra2)
+			for (i in 0:(z$int - 1))
+			{
+				z$sdf[[z$nit + i]] <- z$sdf0[i + 1, , ]
+				if (z$byWave)
+				{
+					z$sdf2[[z$nit + i]] <- z$sdf02[i + 1, , , ]
+				}
+			}
+		}
+		else
+		{
+			if (int==1)
+			{
+				if (!is.null(zz[['sc']]))
+				{
+					z$ssc[z$nit , ,] <- zz$sc
+				}
+			}
+			else
+			{
+                for (i in 1:int)
+                {
+                    if (!is.null(zz[[i]][['sc']]))
+					{
+                        z$ssc[z$nit + (i - 1), , ] <- zz[[i]]$sc
+					}
+                }
+            }
+		}
+		if ((!x$maxlike) && z$cconditional && z$Phase == 3)
+		{
+			if (int == 1)
+			{
+				z$ntim[z$nit, ] <- zz$ntim0
+			}
+			else
+			{
+				for (i in 1:int)
+				{
+					z$ntim[z$nit + (i-1), ] <- zz[[i]]$ntim0
+				}
+			}
+		}
+		## post processing also differs
+		if (phase == 1)
+		{
+			CheckBreaks()
+			if (UserInterruptFlag() || UserRestartFlag())
+			{
+				return(z)
+			}
+			val <- getProgressBar(z$pb)
+			if (z$FinDiff.method)
+			{
+				val <- val + z$pp + 1
+			}
+			else
+			{
+				val <- val + 1
+			}
+			z$pb <- setProgressBar(z$pb, val)
+			progress <- val / z$pb$pbmax * 100
+			if (z$nit <= 5 || z$nit %% z$writefreq == 0 || z$nit %%5 == 0 ||
+				x$maxlike || z$FinDiff.method ||
+				(z$int > 1 && z$nit %% z$writefreq < z$int))
+			{
+				##  Report(c('Phase', z$Phase, 'Iteration ', z$nit, '\n'))
+				if (is.batch())
+				{
+					Report(c('Phase ', z$Phase, ' Iteration ', z$nit,
+							 ' Progress: ', round(progress), '%\n'), sep='')
+				}
+				else
+				{
+					DisplayTheta(z)
+					DisplayDeviations(z, z$sf[z$nit, ])
+				}
+			}
+			if (!z$OK || UserRestartFlag() || UserInterruptFlag())
+				return(z)
+		}
+		else ## phase 3
+		{
+			if (!is.batch())
+			{
+				if (nit < 10)
+				{
+					Report(c("  ", nit, " ",
+							 format(z$sf[nit,], width=1, digits=4, nsmall=4),
+							 "\n"), cf)
+					if (nit >= 10)
+					{
+						CheckBreaks()
+						##    if (nit==10) set up stopkey hint
+						##   Early termination of estimation
+						if( UserInterruptFlag())
+						{
+							Report(c("The user asked for an early stop of the ",
+									 "algorithm during phase 3.\n"), outf)
+							z$Phase3Interrupt <- TRUE
+							if (nit < 500)
+							{
+								if (EarlyEndPhase2Flag())
+								{
+									Report('This implies that ', outf)
+								}
+								else
+								{
+									Report(c("This implies that the estimates",
+											 "are as usual,\nbut the "), outf)
+								}
+								Report(c("diagnostic checks, covariance",
+										 "matrices and t-values \nare less ",
+										 "reliable, because they ",
+										 'are now based on only', nit,
+										 'phase-3 iterations.\n\n'), outf)
+							}
+							z$sf <- z$sf[1:nit, , drop=FALSE]
+							z$sf2 <- z$sf2[1:nit, , , drop=FALSE]
+							if (!x$maxlike)
+							{
+								z$ssc <- z$ssc[1:nit, , , drop=FALSE]
+							}
+							else
+							{
+								z$sdf <-z$sdf[1:nit, , , drop=FALSE]
+								z$sdf2 <-z$sdf2[1:nit, , , ,drop=FALSE]
+							}
+							z$sims <-z$sims[1:nit]
+                            z$Phase3nits <- nit
+                            z$n3 <- nit
+							break
+						}
+						if (UserRestartFlag())
+							break
+					}
+				}
+			}
+			if (!z$OK || UserRestartFlag())
+			{
+				return(z)
+			}
+		}
+	}
+	z
+}
+

Modified: pkg/RSiena/changeLog
===================================================================
--- pkg/RSiena/changeLog	2011-10-28 12:43:23 UTC (rev 179)
+++ pkg/RSiena/changeLog	2011-11-04 12:33:30 UTC (rev 180)
@@ -1,3 +1,11 @@
+2011-11-04 R-forge revision 180.
+
+	* R/phase1.r, R/phase3.r: revision to avoid time consuming copies
+	* src/model/StatisticalCalculator.cpp: stop setting up of networks
+	if no endowment or creation effects.
+	* src/model/effects/behaviorInteractionEffect.cpp: correct properties to be 
+	eol-style not executable
+
 2011-10-28 R-forge revision 179.
 
 	* R/robmon.r: fix problem with forking clusters.

Modified: pkg/RSiena/man/RSiena-package.Rd
===================================================================
--- pkg/RSiena/man/RSiena-package.Rd	2011-10-28 12:43:23 UTC (rev 179)
+++ pkg/RSiena/man/RSiena-package.Rd	2011-11-04 12:33:30 UTC (rev 180)
@@ -30,8 +30,8 @@
 \tabular{ll}{
 Package: \tab RSiena\cr
 Type: \tab Package\cr
-Version: \tab 1.0.12.179\cr
-Date: \tab 2011-10-28\cr
+Version: \tab 1.0.12.180\cr
+Date: \tab 2011-11-04\cr
 License: \tab GPL-2 \cr
 LazyLoad: \tab yes\cr
 }

Modified: pkg/RSiena/src/model/StatisticCalculator.cpp
===================================================================
--- pkg/RSiena/src/model/StatisticCalculator.cpp	2011-10-28 12:43:23 UTC (rev 179)
+++ pkg/RSiena/src/model/StatisticCalculator.cpp	2011-11-04 12:33:30 UTC (rev 180)
@@ -35,6 +35,7 @@
 #include "model/variables/NetworkVariable.h"
 #include "model/variables/BehaviorVariable.h"
 #include "model/tables/Cache.h"
+#include "network/IncidentTieIterator.h"
 
 namespace siena
 {
@@ -293,82 +294,90 @@
 void StatisticCalculator::calculateNetworkEndowmentStatistics(
 	NetworkLongitudinalData * pNetworkData)
 {
-	// In order to calculate the statistics of network endowment effects,
-	// we need the initial network of the given period, and the network
-	// of lost ties, namely, the ties that are present in the initial
-	// network, not missing at the start of the period, and absent in the
-	// current network.
+	// First make sure we have some endowment effects selected
 
-	const Network * pInitialNetwork = pNetworkData->pNetwork(this->lperiod);
-	Network * pLostTieNetwork = pInitialNetwork->clone();
+	const vector<EffectInfo *> & rEffects =
+		this->lpModel->rEndowmentEffects(pNetworkData->name());
 
-	// Duplicate the current network so can overwrite the structurals
-	Network * pCurrentNetwork =
-		this->lpState->pNetwork(pNetworkData->name())->clone();
+    if (rEffects.size() > 0)
+	{
+		// In order to calculate the statistics of network endowment effects,
+		// we need the initial network of the given period, and the network
+		// of lost ties, namely, the ties that are present in the initial
+		// network, not missing at the start of the period, and absent in the
+		// current network.
 
-	// Replace values for structurally determined values from previous
-	// or current period
-	this->replaceNetwork(pCurrentNetwork,
-		pNetworkData->pNetwork(this->lperiod + 1),
-		pNetworkData->pStructuralTieNetwork(this->lperiod + 1));
+		const Network * pInitialNetwork = pNetworkData->pNetwork(this->lperiod);
+		Network * pLostTieNetwork = pInitialNetwork->clone();
 
-	this->replaceNetwork(pCurrentNetwork,
-		pNetworkData->pNetwork(this->lperiod),
-		pNetworkData->pStructuralTieNetwork(this->lperiod));
+		// Duplicate the current network so can overwrite the structurals
+		Network * pCurrentNetwork =
+			this->lpState->pNetwork(pNetworkData->name())->clone();
 
-	// remove missings and current
+		// Replace values for structurally determined values from previous
+		// or current period
+		this->replaceNetwork(pCurrentNetwork,
+			pNetworkData->pNetwork(this->lperiod + 1),
+			pNetworkData->pStructuralTieNetwork(this->lperiod + 1));
 
-	this->subtractNetwork(pLostTieNetwork,
-		pCurrentNetwork);
-	this->subtractNetwork(pLostTieNetwork,
-		pNetworkData->pMissingTieNetwork(this->lperiod + 1));
+		this->replaceNetwork(pCurrentNetwork,
+			pNetworkData->pNetwork(this->lperiod),
+			pNetworkData->pStructuralTieNetwork(this->lperiod));
 
-	// create a new predictor network with only some missings removed
-	Network * pPredictor =
-		pNetworkData->pNetwork(this->lperiod)->clone();
-	this->subtractNetwork(pPredictor,
-		pNetworkData->pMissingTieNetwork(this->lperiod));
+		// remove missings and current
 
-	// We want to pass all networks to the effects in a single state,
-	// hence we overwrite the network in the predictor state.
+		this->subtractNetwork(pLostTieNetwork,
+			pCurrentNetwork);
+		this->subtractNetwork(pLostTieNetwork,
+			pNetworkData->pMissingTieNetwork(this->lperiod + 1));
 
-	string name = pNetworkData->name();
-	const Network * pPredictorNetwork = this->lpPredictorState->pNetwork(name);
-	this->lpPredictorState->pNetwork(name, pPredictor);
+		// create a new predictor network with only some missings removed
+		Network * pPredictor =
+			pNetworkData->pNetwork(this->lperiod)->clone();
+		this->subtractNetwork(pPredictor,
+			pNetworkData->pMissingTieNetwork(this->lperiod));
 
-	// Loop through the endowment effects, calculate the statistics,
-	// and store them.
+		// We want to pass all networks to the effects in a single state,
+		// hence we overwrite the network in the predictor state.
 
-	const vector<EffectInfo *> & rEffects =
-		this->lpModel->rEndowmentEffects(pNetworkData->name());
+		string name = pNetworkData->name();
+		const Network * pPredictorNetwork = this->lpPredictorState->pNetwork(name);
+		this->lpPredictorState->pNetwork(name, pPredictor);
 
- 	EffectFactory factory(this->lpData);
- 	Cache cache;
+		// Loop through the endowment effects, calculate the statistics,
+		// and store them.
 
-	for (unsigned i = 0; i < rEffects.size(); i++)
-	{
-		EffectInfo * pInfo = rEffects[i];
-		NetworkEffect * pEffect =
-			(NetworkEffect *) factory.createEffect(pInfo);
+		//const vector<EffectInfo *> & rEffects =
+		//this->lpModel->rEndowmentEffects(pNetworkData->name());
 
-		// Initialize the effect to work with our data and state of variables.
+		EffectFactory factory(this->lpData);
+		Cache cache;
 
-		pEffect->initialize(this->lpData,
-			this->lpPredictorState,
-			this->lperiod,
-			&cache);
+		for (unsigned i = 0; i < rEffects.size(); i++)
+		{
+			EffectInfo * pInfo = rEffects[i];
+			NetworkEffect * pEffect =
+				(NetworkEffect *) factory.createEffect(pInfo);
 
-		this->lstatistics[pInfo] =
-			pEffect->endowmentStatistic(pLostTieNetwork);
-		delete pEffect;
-	}
+			// Initialize the effect to work with our data and state of variables.
 
-	// Restore the predictor network
-	this->lpPredictorState->pNetwork(name, pPredictorNetwork);
+			pEffect->initialize(this->lpData,
+				this->lpPredictorState,
+				this->lperiod,
+				&cache);
 
-	delete pPredictor;
-	delete pCurrentNetwork;
-	delete pLostTieNetwork;
+			this->lstatistics[pInfo] =
+				pEffect->endowmentStatistic(pLostTieNetwork);
+			delete pEffect;
+		}
+
+		// Restore the predictor network
+		this->lpPredictorState->pNetwork(name, pPredictorNetwork);
+
+		delete pPredictor;
+		delete pCurrentNetwork;
+		delete pLostTieNetwork;
+	}
 }
 
 
@@ -379,6 +388,10 @@
 void StatisticCalculator::calculateNetworkCreationStatistics(
 	NetworkLongitudinalData * pNetworkData)
 {
+ 	const vector<EffectInfo *> & rEffects =
+ 		this->lpModel->rCreationEffects(pNetworkData->name());
+	if (rEffects.size() > 0)
+	{
 	// Duplicate the current network and remove those ties that are
 	// missing at either end of the period.
 
@@ -425,8 +438,8 @@
 	// Loop through the creation effects, calculate the statistics,
 	// and store them.
 
- 	const vector<EffectInfo *> & rEffects =
- 		this->lpModel->rCreationEffects(pNetworkData->name());
+	//	const vector<EffectInfo *> & rEffects =
+ 	//	this->lpModel->rCreationEffects(pNetworkData->name());
 
  	EffectFactory factory(this->lpData);
  	Cache cache;
@@ -455,6 +468,7 @@
 
 	delete pNetwork;
 	delete pGainedTieNetwork;
+	}
 }
 
 
@@ -959,7 +973,7 @@
 				}
 				//}
 		}
-		else
+		else if (rateType == "structural")
 		{
 			// We expect a structural (network-dependent) rate effect here.
 
@@ -1006,6 +1020,48 @@
 			this->lstatistics[pInfo] = statistic;
 			delete pStructural;
 		}
+		else if (rateType == "diffusion")
+		{
+			NetworkLongitudinalData *pNetworkData = this->lpData->
+				pNetworkData(interactionName);
+			Network * pStructural =
+				pNetworkData->pNetwork(this->lperiod)->clone();
+			this->subtractNetwork(pStructural,
+				pNetworkData->pMissingTieNetwork(this->lperiod));
+
+			double statistic = 0;
+			for (int i = 0; i < pBehaviorData->n(); i++)
+			{
+				if (effectName == "avExposure")
+				{
+					double totalAlterValue = 0;
+					double averageAlterValue = 0;
+					if (pStructural->outDegree(i) > 0)
+				    {
+						for (IncidentTieIterator iter = pStructural->outTies(i);
+							 iter.valid();
+							 iter.next())
+						{
+							double alterValue = pBehaviorData->
+								value(this->lperiod,iter.actor());
+							totalAlterValue += alterValue;
+						}
+						averageAlterValue = totalAlterValue /
+							pStructural->outDegree(i);
+				    }
+
+					statistic += averageAlterValue *
+						difference[i];
+				}
+				else
+				{
+				    throw domain_error("Unexpected rate effect " + effectName);
+				}
+			}
+
+			this->lstatistics[pInfo] = statistic;
+			delete pStructural;
+		}
 	}
 
 	delete[] difference;

Modified: pkg/RSiena/src/model/effects/BehaviorInteractionEffect.cpp
===================================================================
--- pkg/RSiena/src/model/effects/BehaviorInteractionEffect.cpp	2011-10-28 12:43:23 UTC (rev 179)
+++ pkg/RSiena/src/model/effects/BehaviorInteractionEffect.cpp	2011-11-04 12:33:30 UTC (rev 180)
@@ -1,177 +1,177 @@
-/******************************************************************************
- * SIENA: Simulation Investigation for Empirical Network Analysis
- *
- * Web: http://www.stats.ox.ac.uk/~snijders/siena/
- *
- * File: BehaviorInteractionEffect.cpp
- *
- * Description: This file contains the implementation of the class
- * BehaviorInteractionEffect.
- *****************************************************************************/
-
-#include "BehaviorInteractionEffect.h"
-
-namespace siena
-{
-
-/**
- * Constructs a new interaction effect between the given effects.
- * The parameter pEffect3 should be 0 for two-way interactions.
- * This effect takes the ownership of the given effects, which mean
- * that the given effects will be destroyed as soon as this
- * effect is destroyed.
- */
-BehaviorInteractionEffect::BehaviorInteractionEffect(
-	const EffectInfo * pEffectInfo,
-	BehaviorEffect * pEffect1,
-	BehaviorEffect * pEffect2,
-	BehaviorEffect * pEffect3) : BehaviorEffect(pEffectInfo)
-{
-	this->lpEffect1 = pEffect1;
-	this->lpEffect2 = pEffect2;
-	this->lpEffect3 = pEffect3;
-}
-
-
-/**
- * Deallocates this effects.
- */
-BehaviorInteractionEffect::~BehaviorInteractionEffect()
-{
-	delete this->lpEffect1;
-	delete this->lpEffect2;
-	delete this->lpEffect3;
-}
-
-
-/**
- * Initializes this effect.
- * @param[in] pData the observed data
- * @param[in] pState the current state of the dependent variables
- * @param[in] period the period of interest
- * @param[in] pCache the cache object to be used to speed up calculations
- */
-void BehaviorInteractionEffect::initialize(const Data * pData,
-	State * pState,
-	int period,
-	Cache * pCache)
-{
-	BehaviorEffect::initialize(pData, pState, period, pCache);
-	this->lpEffect1->initialize(pData, pState, period, pCache);
-	this->lpEffect2->initialize(pData, pState, period, pCache);
-
-	if (this->lpEffect3)
-	{
-		this->lpEffect3->initialize(pData, pState, period, pCache);
-	}
-}
-
-
-/**
- * Does the necessary preprocessing work for calculating the tie flip
- * contributions for a specific ego. This method must be invoked before
- * calling BehaviorEffect::calculateChangeContribution(...).
- */
-void BehaviorInteractionEffect::preprocessEgo(int ego)
-{
-	BehaviorEffect::preprocessEgo(ego);
-
-	this->lpEffect1->preprocessEgo(ego);
-	this->lpEffect2->preprocessEgo(ego);
-
-	if (this->lpEffect3)
-	{
-		this->lpEffect3->preprocessEgo(ego);
-	}
-}
-
-
-/**
- * Calculates the change in the statistic corresponding to this effecf if
- * the given actor were to change his behavior by the given amount.
- */
-double BehaviorInteractionEffect::calculateChangeContribution(int actor,
-	int difference)
-{
-	double contribution =
-		this->lpEffect1->calculateChangeContribution(actor, difference) *
-		this->lpEffect2->calculateChangeContribution(actor, difference);
-
-	contribution /= difference;
-
-	if (this->lpEffect3)
-	{
-		contribution *= this->lpEffect3->calculateChangeContribution(actor,
-			difference) / difference;
-	}
-
-	return contribution;
-}
-
-
-/**
- * Calculates the statistic corresponding to the given ego with respect to the
- * given values of the behavior variable.
- */
-double BehaviorInteractionEffect::egoStatistic(int i,
-	double * currentValues)
-{
-	double statistic;
-
-	if (currentValues[i] != 0)
-	{
[TRUNCATED]

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


More information about the Rsiena-commits mailing list