[Rsiena-commits] r265 - in pkg/RSiena: R src src/model src/model/effects src/model/ml src/model/variables

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Mar 14 13:49:49 CET 2014


Author: natalie
Date: 2014-03-14 13:49:49 +0100 (Fri, 14 Mar 2014)
New Revision: 265

Added:
   pkg/RSiena/R/.Rhistory
   pkg/RSiena/R/getActorStatistics.r
Modified:
   pkg/RSiena/R/initializeFRAN.r
   pkg/RSiena/R/simstatsc.r
   pkg/RSiena/src/model/Model.cpp
   pkg/RSiena/src/model/Model.h
   pkg/RSiena/src/model/StatisticCalculator.cpp
   pkg/RSiena/src/model/StatisticCalculator.h
   pkg/RSiena/src/model/effects/BehaviorEffect.cpp
   pkg/RSiena/src/model/effects/BehaviorEffect.h
   pkg/RSiena/src/model/effects/NetworkEffect.cpp
   pkg/RSiena/src/model/effects/NetworkEffect.h
   pkg/RSiena/src/model/ml/MiniStep.cpp
   pkg/RSiena/src/model/ml/MiniStep.h
   pkg/RSiena/src/model/variables/BehaviorVariable.cpp
   pkg/RSiena/src/model/variables/DependentVariable.cpp
   pkg/RSiena/src/model/variables/DependentVariable.h
   pkg/RSiena/src/model/variables/NetworkVariable.cpp
   pkg/RSiena/src/siena07internals.cpp
   pkg/RSiena/src/siena07internals.h
   pkg/RSiena/src/siena07models.cpp
   pkg/RSiena/src/siena07setup.cpp
   pkg/RSiena/src/siena07utilities.cpp
   pkg/RSiena/src/siena07utilities.h
Log:
Include same changes as in the previous commit in RSienaTest also in RSiena

Added: pkg/RSiena/R/.Rhistory
===================================================================
--- pkg/RSiena/R/.Rhistory	                        (rev 0)
+++ pkg/RSiena/R/.Rhistory	2014-03-14 12:49:49 UTC (rev 265)
@@ -0,0 +1 @@
+print("Hi Natalie")


Property changes on: pkg/RSiena/R/.Rhistory
___________________________________________________________________
Added: svn:mime-type
   + text/plain

Added: pkg/RSiena/R/getActorStatistics.r
===================================================================
--- pkg/RSiena/R/getActorStatistics.r	                        (rev 0)
+++ pkg/RSiena/R/getActorStatistics.r	2014-03-14 12:49:49 UTC (rev 265)
@@ -0,0 +1,59 @@
+##@getActorStatistics. Use as RSiena:::getActorStatistics
+getActorStatistics <- function(algorithm, data, effects)
+{
+	## The following initializations data, effects, and model
+	## for calling "getTargets" in "siena07.setup.h"
+	## is more or less copied from "getTargets" in "getTargets.r".
+	## However, some modifications have been necessary to get it to work.
+	f <- unpackData(data,algorithm)
+	
+	effects <- effects[effects$include,]
+	if (!is.null(algorithm$settings))
+	{
+		effects <- addSettingsEffects(effects, algorithm)
+	}
+	else
+	{
+		effects$setting <- rep("", nrow(effects))
+	}
+	pData <- .Call('setupData', PACKAGE=pkgname,
+			list(as.integer(f$observations)),
+			list(f$nodeSets))
+	## register a finalizer
+	ans <- reg.finalizer(pData, clearData, onexit = FALSE)
+	ans<- .Call('OneMode', PACKAGE=pkgname,
+			pData, list(f$nets))
+	ans<- .Call('Behavior', PACKAGE=pkgname, pData,
+			list(f$behavs))
+	ans<-.Call('ConstantCovariates', PACKAGE=pkgname,
+			pData, list(f$cCovars))
+	ans<-.Call('ChangingCovariates',PACKAGE=pkgname,
+			pData,list(f$vCovars))
+	ans<-.Call('DyadicCovariates',PACKAGE=pkgname,
+			pData,list(f$dycCovars))
+	ans<-.Call('ChangingDyadicCovariates',PACKAGE=pkgname,
+			pData, list(f$dyvCovars))
+	
+	storage.mode(effects$parm) <- 'integer'
+	storage.mode(effects$group) <- 'integer'
+	storage.mode(effects$period) <- 'integer'
+	
+	effects$effectPtr <- rep(NA, nrow(effects))
+	depvarnames <- names(data$depvars)
+	tmpeffects <- split(effects, effects$name)
+	myeffectsOrder <- match(depvarnames, names(tmpeffects))
+	ans <- .Call("effects", PACKAGE=pkgname, pData, tmpeffects)
+	pModel <- ans[[1]][[1]]
+	for (i in 1:length(ans[[2]])) 
+	{
+		effectPtr <- ans[[2]][[i]]
+		tmpeffects[[i]]$effectPtr <- effectPtr
+	}	
+	myeffects <- tmpeffects
+	for(i in 1:length(myeffectsOrder)){
+		myeffects[[i]]<-tmpeffects[[myeffectsOrder[i]]]
+	}
+	ans <- .Call("getTargets", PACKAGE=pkgname, pData, pModel, myeffects, parallelrun=TRUE, returnActorStatistics=TRUE, returnStaticChangeContributions=FALSE)
+	ans	
+}
+

Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/R/initializeFRAN.r	2014-03-14 12:49:49 UTC (rev 265)
@@ -479,7 +479,7 @@
     if (!initC)
     {
         ans <- .Call("getTargets", PACKAGE=pkgname, pData, pModel, myeffects,
-					 z$parallelTesting)
+					 z$parallelTesting, returnActorStatistics=FALSE, returnStaticChangeContributions=FALSE)
 		##stop("done")
 		## create a grid of periods with group names in case want to
 		## parallelize using this or to access chains easily

Modified: pkg/RSiena/R/simstatsc.r
===================================================================
--- pkg/RSiena/R/simstatsc.r	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/R/simstatsc.r	2014-03-14 12:49:49 UTC (rev 265)
@@ -33,6 +33,14 @@
     {
         returnDeps <- z$returnDeps
     }
+	if (is.null(z$returnActorStatistics))
+	{
+		z$returnActorStatistics <- FALSE
+	}
+	if (is.null(z$returnChangeContributions))
+	{
+		z$returnChangeContributions <- FALSE
+	}
     if (is.null(f$seeds))
     {
         seeds <- NULL
@@ -77,7 +85,8 @@
 				 fromFiniteDiff, f$pModel, f$myeffects, z$theta,
 				 randomseed2, returnDeps, z$FinDiff.method,
 				 !is.null(z$cl) && useStreams, z$addChainToStore,
-				  z$returnChains, returnLoglik)
+				  z$returnChains, returnLoglik, 
+				  z$returnActorStatistics, z$returnChangeContributions)
     if (!fromFiniteDiff)
     {
         if (z$FinDiff.method)
@@ -119,6 +128,22 @@
 	{
 		loglik <- NULL
 	}
+	if(z$returnChangeContributions)
+	{	
+		changeContributions <- ans[[9]]
+	}
+	else
+	{
+		changeContributions <- NULL
+	}
+	if(z$returnActorStatistics)
+	{
+		actorStatistics <- ans[[10]]
+	}
+	else
+	{
+		actorStatistics <- NULL
+	}
     if (returnDeps)
     {
         ## attach the names
@@ -137,7 +162,8 @@
     }
 		 ## browser()
     list(sc = sc, fra = fra, ntim0 = ntim, feasible = TRUE, OK = TRUE,
-         sims=sims, f$seeds, chain=chain, loglik=loglik)
+         sims=sims, f$seeds, chain=chain, loglik=loglik,
+		 actorStatistics = actorStatistics, changeContributions = changeContributions)
 }
 
 ##@clearData siena07 Finalizer to clear Data object in C++

Modified: pkg/RSiena/src/model/Model.cpp
===================================================================
--- pkg/RSiena/src/model/Model.cpp	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/Model.cpp	2014-03-14 12:49:49 UTC (rev 265)
@@ -195,6 +195,22 @@
 }
 
 /**
+ * Stores if change contribution are needed in the current simulation
+ */
+void Model::needChangeContributions(bool flag)
+{
+	this->lneedChangeContributions2 = flag;
+}
+
+/**
+ * Returns if change contribution are needed in the current simulation
+ */
+bool Model::needChangeContributions() const
+{
+	return this->lneedChangeContributions2;
+}
+
+/**
  * Stores the number of ML steps
  */
 void Model::numberMLSteps(int value)

Modified: pkg/RSiena/src/model/Model.h
===================================================================
--- pkg/RSiena/src/model/Model.h	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/Model.h	2014-03-14 12:49:49 UTC (rev 265)
@@ -130,6 +130,9 @@
 	void parallelRun(bool flag);
 	bool parallelRun() const;
 
+	void needChangeContributions(bool flag);
+	bool needChangeContributions() const;
+
 	void modelType(int type);
 	ModelType modelType() const;
 	bool modelTypeB() const;
@@ -244,6 +247,9 @@
 	// and score calculations
 	bool lparallelRun;
 
+	//indicates whether change contributions are needed
+	bool lneedChangeContributions2;
+
 	// number of steps in a run for ML
 	int lnumberMLSteps;
 

Modified: pkg/RSiena/src/model/StatisticCalculator.cpp
===================================================================
--- pkg/RSiena/src/model/StatisticCalculator.cpp	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/StatisticCalculator.cpp	2014-03-14 12:49:49 UTC (rev 265)
@@ -13,6 +13,8 @@
 #include <stdexcept>
 #include <vector>
 #include <cmath>
+#include <R_ext/Arith.h>
+#include <R_ext/Print.h>
 
 #include "StatisticCalculator.h"
 #include "data/Data.h"
@@ -65,12 +67,59 @@
 	this->lperiod = period;
 	this->lpPredictorState = new State();
 	this->lpStateLessMissingsEtc = new State();
+	this->lneedActorStatistics = 0;
+	this->lcountStaticChangeContributions = 0;
 
 	this->calculateStatistics();
 }
 
+/**
+ * Constructor.
+ * @param[in] pData the observed data
+ * @param[in] pModel the model whose effect statistics are to be calculated
+ * @param[in] pState the current state of the dependent variables
+ * @param[in] period the period under consideration
+ * @param[in] returnActorStatistics whether individual actor statistics should be returned
+ */
+StatisticCalculator::StatisticCalculator(const Data * pData, const Model * pModel, State * pState, int period, bool returnActorStatistics)
+{
+	this->lpData = pData;
+	this->lpModel = pModel;
+	this->lpState = pState;
+	this->lperiod = period;
+	this->lpPredictorState = new State();
+	this->lpStateLessMissingsEtc = new State();
+	this->lneedActorStatistics = returnActorStatistics;
+	this->lcountStaticChangeContributions = 0;
 
+	this->calculateStatistics();
+}
+
 /**
+ * Constructor.
+ * @param[in] pData the observed data
+ * @param[in] pModel the model whose effect statistics are to be calculated
+ * @param[in] pState the current state of the dependent variables
+ * @param[in] period the period under consideration
+ * @param[in] returnActorStatistics whether individual actor statistics should be returned
+ * @param[in] returnStaticChangeContributions whether contributions of effects on possible next tie flips or behavior changes are needed
+ */
+StatisticCalculator::StatisticCalculator(const Data * pData, const Model * pModel, State * pState, int period, bool returnActorStatistics, bool returnStaticChangeContributions)
+{
+	this->lpData = pData;
+	this->lpModel = pModel;
+	this->lpState = pState;
+	this->lperiod = period;
+	this->lpPredictorState = new State();
+	this->lpStateLessMissingsEtc = new State();
+	this->lneedActorStatistics = returnActorStatistics;
+	this->lcountStaticChangeContributions = returnStaticChangeContributions;
+
+	this->calculateStatistics();
+}
+
+
+/**
  * Deallocates this model.
  */
 StatisticCalculator::~StatisticCalculator()
@@ -109,6 +158,24 @@
 	this->lpStateLessMissingsEtc->deleteValues();
 	delete this->lpStateLessMissingsEtc;
 	this->lpStateLessMissingsEtc = 0;
+
+	while (!this->lstaticChangeContributions.empty())
+	{
+		vector<double *> vec = this->lstaticChangeContributions.begin()->second;
+		while (!vec.empty())
+		{
+			double * array = *(vec.begin());
+			vec.erase(vec.begin());
+			delete array;
+		}
+		this->lstaticChangeContributions.erase(this->lstaticChangeContributions.begin());
+	}
+	while (!this->lactorStatistics.empty())
+	{
+		double * array = this->lactorStatistics.begin()->second;
+		this->lactorStatistics.erase(this->lactorStatistics.begin());
+		delete[] array;
+	}
 }
 
 
@@ -133,8 +200,38 @@
 	return iter->second;
 }
 
+/**
+ * Returns the tie flip contributions or the behavior change contributions of the given effect.
+ */
+vector<double *> StatisticCalculator::staticChangeContributions(EffectInfo * pEffect) const
+{
+	map<EffectInfo *, vector<double *> >::const_iterator iter =
+		this->lstaticChangeContributions.find(pEffect);
+	if (iter == this->lstaticChangeContributions.end())
+	{
+		throw invalid_argument(
+			"Unknown effect: The given effect is not part of the model.");
+	}
+	return iter->second;
+}
 
 /**
+ * Returns the actor statistics of the given effect.
+ */
+double * StatisticCalculator::actorStatistics(EffectInfo * pEffect) const
+{
+	map<EffectInfo *, double * >::const_iterator iter =
+		this->lactorStatistics.find(pEffect);
+
+	if (iter == this->lactorStatistics.end())
+	{
+		throw invalid_argument(
+			"Unknown effect: The given effect is not part of the model.");
+	}
+	return iter->second;
+}
+
+/**
  * Returns the simulated distance for the given network and period.
  */
 int StatisticCalculator::distance(LongitudinalData * pData, int period)
@@ -312,7 +409,51 @@
 			this->lperiod,
 			&cache);
 
-		this->lstatistics[pInfo] = pEffect->evaluationStatistic();
+		if(this->lneedActorStatistics)
+		{
+			pair<double, double * > p = pEffect->evaluationStatistic(this->lneedActorStatistics);
+			this->lstatistics[pInfo] = p.first;
+			this->lactorStatistics[pInfo] = p.second;
+		}
+		else
+		{
+			this->lstatistics[pInfo] = pEffect->evaluationStatistic();
+		}
+		if(this->lcountStaticChangeContributions)
+		{
+			int egos =  pCurrentLessMissingsEtc->n();
+			int alters =  pCurrentLessMissingsEtc->m();
+			vector<double *> egosMap(egos);
+			this->lstaticChangeContributions.insert(make_pair(pInfo,egosMap));
+			for (int e = 0; e < egos ; e++)
+			{
+				cache.initialize(e);
+				pEffect->initialize(this->lpData, this->lpPredictorState, this->lperiod, &cache);
+				double * contributions = new double[egos];
+				this->lstaticChangeContributions.at(pInfo).at(e) = contributions;
+				pEffect->preprocessEgo(e);
+				//TODO determine permissible changes (see: NetworkVariable::calculatePermissibleChanges())
+				for (int a = 0; a < alters ; a++)
+				{
+					if(a == e)
+					{
+						this->lstaticChangeContributions.at(pInfo).at(e)[a] = 0;
+					}
+					else
+					{
+						// Tie withdrawal contributes the opposite of tie creating
+						if (pCurrentLessMissingsEtc->tieValue(e,a))
+						{
+							this->lstaticChangeContributions.at(pInfo).at(e)[a] = -pEffect->calculateContribution(a);
+						}
+						else
+						{
+							this->lstaticChangeContributions.at(pInfo).at(e)[a] = pEffect->calculateContribution(a);
+						}
+					}
+				}
+			}
+		}
 		delete pEffect;
 	}
 
@@ -397,8 +538,17 @@
 				this->lperiod,
 				&cache);
 
-			this->lstatistics[pInfo] =
-				pEffect->endowmentStatistic(pLostTieNetwork);
+			if(this->lneedActorStatistics)
+			{
+				pair<double, double * > p = pEffect->endowmentStatistic(pLostTieNetwork, this->lneedActorStatistics);
+				this->lstatistics[pInfo] = p.first;
+				this->lactorStatistics[pInfo] = p.second;
+			}
+			else
+			{
+				this->lstatistics[pInfo] =
+					pEffect->endowmentStatistic(pLostTieNetwork);
+			}
 			delete pEffect;
 		}
 
@@ -467,9 +617,17 @@
 				this->lperiod,
 				&cache);
 
-			this->lstatistics[pInfo] =
-				pEffect->creationStatistic(pGainedTieNetwork);
-
+			if(this->lneedActorStatistics)
+			{
+				pair<double, double * > p = pEffect->creationStatistic(pGainedTieNetwork, this->lneedActorStatistics);
+				this->lstatistics[pInfo] = p.first;
+				this->lactorStatistics[pInfo] = p.second;
+			}
+			else
+			{
+				this->lstatistics[pInfo] =
+					pEffect->creationStatistic(pGainedTieNetwork);
+			}
 			delete pEffect;
 		}
 
@@ -544,8 +702,51 @@
 			this->lperiod,
 			&cache);
 
-		this->lstatistics[pInfo] =
-			pEffect->evaluationStatistic(currentValues);
+		if(this->lneedActorStatistics)
+		{
+			pair<double, double * > p = pEffect->evaluationStatistic(currentValues,this->lneedActorStatistics);
+			this->lstatistics[pInfo] = p.first;
+			this->lactorStatistics[pInfo] = p.second;
+		}
+		else
+		{
+			this->lstatistics[pInfo] =
+				pEffect->evaluationStatistic(currentValues);
+		}
+		if(this->lcountStaticChangeContributions)
+		{
+			int  choices = 3;
+			vector<double *> egosMap(pBehaviorData->n());
+			this->lstaticChangeContributions.insert(make_pair(pInfo,egosMap));
+			for (int e = 0; e < pBehaviorData->n(); e++)
+			{
+				cache.initialize(e);
+			    pEffect->initialize(this->lpData, this->lpPredictorState, this->lperiod, &cache);
+				double * contributions = new double[choices];
+				this->lstaticChangeContributions.at(pInfo).at(e) = contributions;
+				pEffect->preprocessEgo(e);
+				// no change gives no contribution
+				this->lstaticChangeContributions.at(pInfo).at(e)[1] = 0;
+				// calculate the contribution of downward change
+				if (currentValues[e] > pBehaviorData->min() && !pBehaviorData->upOnly(this->lperiod))
+				{
+					this->lstaticChangeContributions.at(pInfo).at(e)[0] = pEffect->calculateChangeContribution(e,-1);
+				}
+				else
+				{
+					this->lstaticChangeContributions.at(pInfo).at(e)[0] = R_NaN;
+				}
+				// calculate the contribution of upward change
+				if (currentValues[e] < pBehaviorData->max() && !pBehaviorData->downOnly(this->lperiod))
+				{
+					this->lstaticChangeContributions.at(pInfo).at(e)[2] = pEffect->calculateChangeContribution(e,1);
+				}
+				else
+				{
+					this->lstaticChangeContributions.at(pInfo).at(e)[2] = R_NaN;
+				}
+			}
+		}
 
 		delete pEffect;
 	}
@@ -569,9 +770,17 @@
 			this->lperiod,
 			&cache);
 
-		this->lstatistics[pInfo] =
-			pEffect->endowmentStatistic(difference, currentValues);
-
+		if(this->lneedActorStatistics)
+		{
+			pair<double, double * > p = pEffect->endowmentStatistic(difference, currentValues,this->lneedActorStatistics);
+			this->lstatistics[pInfo] = p.first;
+			this->lactorStatistics[pInfo] = p.second;
+		}
+		else
+		{
+			this->lstatistics[pInfo] =
+				pEffect->endowmentStatistic(difference, currentValues);
+		}
 		delete pEffect;
 	}
 
@@ -593,10 +802,17 @@
 			this->lpPredictorState,
 			this->lperiod,
 			&cache);
-
-		this->lstatistics[pInfo] =
-			pEffect->creationStatistic(difference, currentValues);
-
+		if(this->lneedActorStatistics)
+		{
+			pair<double, double * > p =  pEffect->creationStatistic(difference, currentValues, this->lneedActorStatistics);
+			this->lstatistics[pInfo] = p.first;
+			this->lactorStatistics[pInfo] = p.second;
+		}
+		else
+		{
+			this->lstatistics[pInfo] =
+				pEffect->creationStatistic(difference, currentValues);
+		}
 		delete pEffect;
 	}
 

Modified: pkg/RSiena/src/model/StatisticCalculator.h
===================================================================
--- pkg/RSiena/src/model/StatisticCalculator.h	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/StatisticCalculator.h	2014-03-14 12:49:49 UTC (rev 265)
@@ -58,9 +58,17 @@
 		const Model * pModel,
 		State * pState,
 		int period);
+	StatisticCalculator(const Data * pData,
+		const Model * pModel, State * pState,
+		int period, bool returnActorStatistics);
+	StatisticCalculator(const Data * pData,
+		const Model * pModel, State * pState,
+		int period, bool returnActorStatistics, bool returnStaticChangeContributions);
 	virtual ~StatisticCalculator();
 
 	double statistic(EffectInfo * pEffectInfo) const;
+	vector<double *> staticChangeContributions(EffectInfo * pEffect) const;
+	double * actorStatistics(EffectInfo * pEffect) const;
 	int distance(LongitudinalData * pData, int period) const;
 	int settingDistance(LongitudinalData * pData, string setting,
 		int period) const;
@@ -100,9 +108,21 @@
 	// The period of the evolution we are interested in
 	int lperiod;
 
+	// indicates whether actor statistics are needed
+	bool lneedActorStatistics;
+
+	// indicates whether change contributions are needed
+	bool lcountStaticChangeContributions;
+
 	// The resulting map of statistic values
 	map<EffectInfo *, double> lstatistics;
 
+	// The resulting map of actor statistic values
+	map<EffectInfo *, double * > lactorStatistics;
+
+	// The change contributions of all effects
+	map<EffectInfo *, vector<double *> > lstaticChangeContributions;
+
 	// Array of simulated distances per variable
 	map<LongitudinalData *, int *> ldistances;
 

Modified: pkg/RSiena/src/model/effects/BehaviorEffect.cpp
===================================================================
--- pkg/RSiena/src/model/effects/BehaviorEffect.cpp	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/effects/BehaviorEffect.cpp	2014-03-14 12:49:49 UTC (rev 265)
@@ -134,20 +134,38 @@
  */
 double BehaviorEffect::evaluationStatistic(double * currentValues)
 {
+	return this->evaluationStatistic(currentValues, false).first;
+}
+
+pair<double,double *> BehaviorEffect::evaluationStatistic(double * currentValues, bool needActorStatistics)
+{
 	double statistic = 0;
 	int n = this->n();
 
+	double * actorStatistics = 0;
+	if(needActorStatistics)
+	{
+		 actorStatistics = new double[n];
+	}
+
 	for (int i = 0; i < n; i++)
 	{
 		this->preprocessEgo(i);
 		if (!this->missing(this->period(), i) &&
 			!this->missing(this->period() + 1, i))
 		{
-			statistic += this->egoStatistic(i, currentValues);
+			if(needActorStatistics)
+			{
+				actorStatistics[i] = this->egoStatistic(i, currentValues);
+				statistic += actorStatistics[i];
+			}
+			else
+			{
+				statistic += this->egoStatistic(i, currentValues);
+			}
 		}
 	}
-
-	return statistic;
+	return make_pair(statistic,actorStatistics);
 }
 
 
@@ -170,20 +188,37 @@
 double BehaviorEffect::endowmentStatistic(const int * difference,
 	double * currentValues)
 {
+	return endowmentStatistic(difference, currentValues, false).first;
+}
+
+pair<double, double * > BehaviorEffect::endowmentStatistic(const int * difference, double * currentValues, bool needActorStatistics)
+{
 	double statistic = 0;
 	int n = this->n();
 
+	double * actorStatistics = 0;
+	if(needActorStatistics)
+	{
+		actorStatistics = new double[n];
+	}
+
 	for (int i = 0; i < n; i++)
 	{
 		this->preprocessEgo(i);
 		if (!this->missing(this->period(), i))
 		{
-			statistic += this->egoEndowmentStatistic(i, difference,
-				currentValues);
+			if(needActorStatistics)
+			{
+				actorStatistics[i] = this->egoEndowmentStatistic(i, difference, currentValues);
+				statistic += actorStatistics[i];
+			}
+			else
+			{
+				statistic += this->egoEndowmentStatistic(i, difference, currentValues);
+			}
 		}
 	}
-
-	return statistic;
+	return make_pair(statistic, actorStatistics);
 }
 /**
  * Returns the statistic corresponding the given ego as part of
@@ -210,6 +245,12 @@
 double BehaviorEffect::creationStatistic(int * difference,
 	double *currentValues)
 {
+	return creationStatistic(difference, currentValues, false).first;
+}
+
+pair<double, double * > BehaviorEffect::creationStatistic(int * difference,
+	double *currentValues, bool needActorStatistics)
+{
 	// Here we use a trick. The creation statistics are very similar to the
 	// endowmnent statistics, but instead of summing over all actors with
 	// decreasing values, we must now sum over all actors with increasing
@@ -217,23 +258,36 @@
 	// statistic.
 
 	int n = this->n();
-
+	double statistic = 0;
+	double * actorStatistics = 0;
 	for (int i = 0; i < n; i++)
 	{
 		difference[i]  = -difference[i];
 	}
 
-	double statistic = this->endowmentStatistic(difference, currentValues);
+	if(needActorStatistics)
+	{
+		pair<double, double *> p = this->endowmentStatistic(difference, currentValues, needActorStatistics);
+		statistic = p.first;
+		actorStatistics = p.second;
+		for (int i = 0; i < n; i++)
+		{
+			actorStatistics[i]  = -actorStatistics[i];
+		}
+	}
+	else
+	{
+		statistic = this->endowmentStatistic(difference, currentValues);
+	}
 
 	for (int i = 0; i < n; i++)
 	{
 		difference[i]  = -difference[i];
 	}
 
-	return -statistic;
+	return make_pair(-statistic, actorStatistics);
 }
 
-
 /**
  * Does the necessary preprocessing work for calculating the probabilities
  * for a specific ego. This method must be invoked before

Modified: pkg/RSiena/src/model/effects/BehaviorEffect.h
===================================================================
--- pkg/RSiena/src/model/effects/BehaviorEffect.h	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/effects/BehaviorEffect.h	2014-03-14 12:49:49 UTC (rev 265)
@@ -13,7 +13,10 @@
 #define BEHAVIOREFFECT_H_
 
 #include "Effect.h"
+#include <utility>
 
+using namespace std;
+
 namespace siena
 {
 
@@ -51,10 +54,15 @@
 		int difference) = 0;
 
 	virtual double evaluationStatistic(double * currentValues);
+	virtual pair<double,double *> evaluationStatistic(double * currentValues, bool needActorStatistics);
 	virtual double endowmentStatistic(const int * difference,
 		double * currentValues);
+	virtual pair<double,double *> endowmentStatistic(const int * difference,
+		double * currentValues, bool needActorStatistics);
 	virtual double creationStatistic(int * difference,
 		double * currentValues);
+	virtual pair<double,double *> creationStatistic(int * difference,
+		double * currentValues, bool needActorStatistics);
 	virtual double egoStatistic(int ego, double * currentValues);
 	virtual double egoEndowmentStatistic(int ego, const int * difference,
 		double * currentValues);

Modified: pkg/RSiena/src/model/effects/NetworkEffect.cpp
===================================================================
--- pkg/RSiena/src/model/effects/NetworkEffect.cpp	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/effects/NetworkEffect.cpp	2014-03-14 12:49:49 UTC (rev 265)
@@ -127,7 +127,16 @@
 	return this->statistic(this->pNetwork());
 }
 
+/**
+ * Returns the statistic corresponding to this effect as part of
+ * the evaluation function.
+ */
+pair<double, double * > NetworkEffect::evaluationStatistic(bool needActorStatistics)
+{
+		return this->statistic(this->pNetwork(), needActorStatistics);
+}
 
+
 /**
  * Returns the statistic corresponding to this effect as part of
  * the endowment function.
@@ -137,7 +146,16 @@
 	return this->statistic(pLostTieNetwork);
 }
 
+/**
+ * Returns the statistic corresponding to this effect as part of
+ * the endowment function.
+ */
+pair<double, double *> NetworkEffect::endowmentStatistic(Network * pLostTieNetwork, bool needActorStatistics)
+{
+	return this->statistic(pLostTieNetwork, needActorStatistics);
+}
 
+
 /**
  * Returns the statistic corresponding to this effect as part of
  * the creation function.
@@ -152,7 +170,17 @@
 	return this->endowmentStatistic(pGainedTieNetwork);
 }
 
+/**
+ * Returns the statistic corresponding to this effect as part of
+ * the creation function.
+ */
+pair<double, double *> NetworkEffect::creationStatistic(Network * pGainedTieNetwork, bool needActorStatistics)
+{
+	return this->endowmentStatistic(pGainedTieNetwork, needActorStatistics);
+}
 
+
+
 /**
  * Returns if this effect is an ego effect.
  */
@@ -175,24 +203,41 @@
  */
 double NetworkEffect::statistic(const Network * pSummationTieNetwork)
 {
+	return statistic(pSummationTieNetwork, false).first;
+}
+
+pair<double,double *> NetworkEffect::statistic(const Network * pSummationTieNetwork, bool needActorStatistics)
+{
 	this->initializeStatisticCalculation();
 
 	int n = pSummationTieNetwork->n();
 	Cache * pCache = this->pCache();
 	double statistic = 0;
-
+	double * actorStatistics = 0;
+	if(needActorStatistics)
+	{
+	  actorStatistics = new double[n];
+	}
 	for (int i = 0; i < n; i++)
 	{
 		pCache->initialize(i);
 		this->preprocessEgo(i);
 		this->onNextEgo(i);
-		statistic += this->egoStatistic(i, pSummationTieNetwork);
+		if(needActorStatistics)
+		{
+			actorStatistics[i] = this->egoStatistic(i, pSummationTieNetwork);
+			statistic += actorStatistics[i];
+		}
+		else
+		{
+			statistic += this->egoStatistic(i, pSummationTieNetwork);
+		}
 	}
 
 	this->cleanupStatisticCalculation();
 	//Rprintf(" %f sum \n ", statistic);
 
-	return statistic;
+	return make_pair(statistic,actorStatistics);
 }
 
 

Modified: pkg/RSiena/src/model/effects/NetworkEffect.h
===================================================================
--- pkg/RSiena/src/model/effects/NetworkEffect.h	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/effects/NetworkEffect.h	2014-03-14 12:49:49 UTC (rev 265)
@@ -13,7 +13,10 @@
 #define NETWORKEFFECT_H_
 
 #include "Effect.h"
+#include <utility>
 
+using namespace std;
+
 namespace siena
 {
 
@@ -63,13 +66,17 @@
 	virtual double calculateContribution(int alter) const = 0;
 
 	virtual double evaluationStatistic();
+	virtual pair <double, double * > evaluationStatistic(bool needActorStatistics);
 	virtual double endowmentStatistic(Network * pLostTieNetwork);
+	virtual pair <double, double * > endowmentStatistic(Network * pLostTieNetwork, bool needActorStatistics);
 	virtual double creationStatistic(Network * pGainedTieNetwork);
+	virtual pair <double, double * > creationStatistic(Network * pGainedTieNetwork, bool needActorStatistics);
 
 	virtual bool egoEffect() const;
 
 protected:
 	virtual double statistic(const Network * pSummationTieNetwork);
+	virtual pair <double, double * > statistic(const Network * pSummationTieNetwork, bool needActorStatistics);
 	virtual void initializeStatisticCalculation();
 	virtual void onNextEgo(int ego);
 	virtual double egoStatistic(int ego,

Modified: pkg/RSiena/src/model/ml/MiniStep.cpp
===================================================================
--- pkg/RSiena/src/model/ml/MiniStep.cpp	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/ml/MiniStep.cpp	2014-03-14 12:49:49 UTC (rev 265)
@@ -41,6 +41,7 @@
 	this->lmissingIndex = -1;
 	this->lorderingKey = 0;
 	this->ldiagonal = false;
+	this->lpChangeContributions = 0;
 }
 
 
@@ -55,6 +56,10 @@
 	}
 
 	this->lpOption = 0;
+	if (this->lpChangeContributions)
+	{
+		delete this->lpChangeContributions;
+	}
 }
 
 
@@ -281,4 +286,22 @@
 		!missing;
 }
 
+/**
+ * Stores the contributions of each effect to possible
+ * tie flips or behavior changes before this ministep.
+ */
+void MiniStep::changeContributions(map<const EffectInfo *, vector<double> > * contributions)
+{
+	this->lpChangeContributions = contributions;
 }
+
+/**
+ * Returns the contributions of each effect to possible
+ * tie flips or behavior changes before this ministep.
+ */
+map<const EffectInfo *, vector<double> > * MiniStep::changeContributions() const
+{
+	return this->lpChangeContributions;
+}
+
+}

Modified: pkg/RSiena/src/model/ml/MiniStep.h
===================================================================
--- pkg/RSiena/src/model/ml/MiniStep.h	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/ml/MiniStep.h	2014-03-14 12:49:49 UTC (rev 265)
@@ -13,6 +13,8 @@
 #define MINISTEP_H_
 
 #include <string>
+#include <vector>
+#include <map>
 
 using namespace std;
 
@@ -27,6 +29,7 @@
 class DependentVariable;
 class Chain;
 class Option;
+class EffectInfo;
 
 
 // ----------------------------------------------------------------------------
@@ -100,6 +103,9 @@
 
 	virtual bool firstOfConsecutiveCancelingPair() const;
 
+	map<const EffectInfo *, vector<double> >* changeContributions() const;
+	void changeContributions(map<const EffectInfo *, vector<double> > * contributions);
+
 protected:
 	void pOption(const Option * pOption);
     void diagonal(bool value);
@@ -172,6 +178,9 @@
 	// of the same chain and x precedes y in the chain.
 
 	double lorderingKey;
+
+	// Stores for each effect its contributions to the tie flip probabilities or behavior change probabilities
+	map<const EffectInfo *, vector<double> > * lpChangeContributions;
 };
 
 

Modified: pkg/RSiena/src/model/variables/BehaviorVariable.cpp
===================================================================
--- pkg/RSiena/src/model/variables/BehaviorVariable.cpp	2014-03-14 12:49:24 UTC (rev 264)
+++ pkg/RSiena/src/model/variables/BehaviorVariable.cpp	2014-03-14 12:49:49 UTC (rev 265)
@@ -90,6 +90,11 @@
 	this->lendowmentEffectContribution = 0;
 	this->lcreationEffectContribution = 0;
 	this->lprobabilities = 0;
+
+	if(this->lpChangeContribution != 0)
+	{
+		delete this->lpChangeContribution;
+	}
 }
 
 
@@ -268,6 +273,10 @@
 		// insert ministep in chain
 		BehaviorChange * pMiniStep =
 			new BehaviorChange(this->lpData, actor, difference);
+		if (this->pSimulation()->pModel()->needChangeContributions())
+		{
+			pMiniStep->changeContributions(lpChangeContribution);
+		}
 		this->pSimulation()->pChain()->insertBefore(pMiniStep,
 			this->pSimulation()->pChain()->pLast());
 		pMiniStep->logChoiceProbability(log(this->lprobabilities[difference
@@ -310,18 +319,22 @@
 	this->lupPossible = true;
 	this->ldownPossible = true;
 
+	int evaluationEffectCount = this->pEvaluationFunction()->rEffects().size();
+	int endowmentEffectCount = this->pEndowmentFunction()->rEffects().size();
+	int creationEffectCount = this->pCreationFunction()->rEffects().size();
+
 	// initialize for later use!
-	for (unsigned i = 0; i < pEvaluationFunction()->rEffects().size(); i++)
+	for (unsigned i = 0; i < evaluationEffectCount; i++)
 	{
 		this->levaluationEffectContribution[1][i] =	0;
 	}
-	for (unsigned i = 0; i < pEndowmentFunction()->rEffects().size(); i++)
+	for (unsigned i = 0; i < endowmentEffectCount; i++)
 	{
 		this->lendowmentEffectContribution[1][i] =	0;
 		this->lendowmentEffectContribution[2][i] = 	0;
 	}
 
-	for (unsigned i = 0; i < this->pCreationFunction()->rEffects().size(); i++)
+	for (unsigned i = 0; i < creationEffectCount; i++)
 	{
 		for (int j = 0; j < 3; j++)
 		{
@@ -329,6 +342,26 @@
 		}
 	}
 
+	if (this->pSimulation()->pModel()->needChangeContributions())
+	{
+		this->lpChangeContribution = new map<const EffectInfo *, vector<double> >[evaluationEffectCount+endowmentEffectCount+creationEffectCount];
+		for (int i = 0; i < evaluationEffectCount; i++)
+		{
+			vector<double> vec(3,0);
+			this->lpChangeContribution->insert(make_pair(this->pEvaluationFunction()->rEffects()[i]->pEffectInfo(), vec));
+		}
+		for (int i = 0; i < endowmentEffectCount; i++)
+		{
+			vector<double> vec(3,0);
+			this->lpChangeContribution->insert(make_pair(this->pEndowmentFunction()->rEffects()[i]->pEffectInfo(), vec));
+		}
+		for (int i = 0; i < creationEffectCount; i++)
+		{
+			vector<double> vec(3,0);
+			this->lpChangeContribution->insert(make_pair(this->pCreationFunction()->rEffects()[i]->pEffectInfo(), vec));
+		}
+	}
+
 	// Calculate the probability for downward change.
 	// Defer exp till can subtract the largest to avoid overflow.
 
@@ -432,6 +465,10 @@
 			(BehaviorEffect *) pFunction->rEffects()[i];
 		double thisContribution =
 			pEffect->calculateChangeContribution(actor, difference);
+		if (this->pSimulation()->pModel()->needChangeContributions())
+		{
+			(* this->lpChangeContribution).at(pEffect->pEffectInfo()).at(difference + 1) = thisContribution;
+		}
 		this->levaluationEffectContribution[difference + 1][i] =
 			thisContribution;
 		contribution += pEffect->parameter() * thisContribution;
@@ -451,6 +488,10 @@
 			(BehaviorEffect *) pFunction->rEffects()[i];
 		double thisContribution =
 			pEffect->calculateChangeContribution(actor, difference);
+		if (this->pSimulation()->pModel()->needChangeContributions())
+		{
+			(* this->lpChangeContribution).at(pEffect->pEffectInfo()).at(difference + 1) = thisContribution;
+		}
 		this->lendowmentEffectContribution[difference + 1][i] =
 			thisContribution;
 		contribution += pEffect->parameter() * thisContribution;
@@ -476,6 +517,10 @@
 			(BehaviorEffect *) pFunction->rEffects()[i];
 		double thisContribution =
 			pEffect->calculateChangeContribution(actor, difference);
+		if (this->pSimulation()->pModel()->needChangeContributions())
+		{
+			(* this->lpChangeContribution).at(pEffect->pEffectInfo()).at(difference + 1) = thisContribution;
+		}
 		this->lcreationEffectContribution[difference + 1][i] =
 			thisContribution;
 		contribution += pEffect->parameter() * thisContribution;

Modified: pkg/RSiena/src/model/variables/DependentVariable.cpp
[TRUNCATED]

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


More information about the Rsiena-commits mailing list