Infected compartment spatRaster layer is severely bugged

[bryce-carson]

Background

During testing of recent documentation changes, and a small readability refactor, a fatal error condition I had written surfaced. During inspection I found that the Infected layer is wildly incorrect, with cells exhibiting a strange latitudinal pattern.

I need to resolve this, and I imagine that once I do it will also resolve issue #24.

Susceptible

Infected

Population

This plot shows that a newly acquired raster of class SpatRaster with no aggregation does not actually contain NAs where I thought it did. The values in the raster outside of what is the Democratic Republic of the Congo (the spatial domain of the DRC), are also zero, like zero values within its spatial domain, not NA.

cod <- getCountryPopulation.SpatRaster("COD")
terra::plot(terra::classify(cod, cuts, right = FALSE), colNA = "black")
any(is.na(values(cod))) # FALSE

[bryce-carson]

Investigation

These two plots of the SVEIRD layers for Ituri and Nord-Kivu display that the areas outside of the spatial domain (outside of the two provinces) are indeed NA, which is what I'd expect given the following code.

  population <- terra::crop(population, subregions, mask = TRUE)

The function documentation for terra::crop states that the mask parameter, which I use, will cause values outside the lay of land within the the masking SpatVector also used to crop will become NA (what I want and expect).

If ‘x’ is a ‘SpatRaster’: Create a new SpatRaster that has the same values as SpatRaster ‘x’, except for the cells that are ‘NA’ (or other ‘maskvalue’) in another SpatRaster (the 'mask'), or the cells that are not covered by a SpatVector or SpatExtent. These cells become ‘NA’ (or another ‘updatevalue’).

Conclusion

I can conclude that the issue wherein values of the Infected layer become strange is due to a function which modifies the layer, not creates it. This leaves the following guideline to find the issue: any function or expression which runs and causes Infected to be "funny" is the culprit. I'm betting it is either assimilate data or the basic process implementing the compartmental model and which overwrites the layers each loop iteration.

[bryce-carson]

The layer looked fine until after the first occurrence of data assimilation, so something is implemented incorrectly in that process which absolutely must be resolved for anything whatsoever to work.

However, data assimilation does not occur in a non-BDA simulation, which is what I opened the other issue with and what I believe I witness the first example of the funny business in the infected compartment layer.

[bryce-carson]

This image shows the extent of the bug(s) as time progresses. Each subplot is a single timestep, with the first image being the Infected layer at the beginning of the simulation, immediately after seeding, which I've determined is quite correct (further testing regarding extent required).

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 924 to 925 in 73268b8

	layers %<>% castSeedDataQueensNeighbourhood(seedData, neighbourhood.order)
	timeseries <- terra::sds("names<-"(as.list(layers), names(layers)))

The following subplots are the Infected layer at the end of each timestep (the second subplot being at the end of the first timestep, actually, since the first subplot is t=0 [prior to the loop]).

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 1062 to 1065 in 0df1abf

	for(layer in names(layers)) {
	i <- which(layer == names(layers))
	terra::add(timeseries[i]) <- layers[layer]
	}

[bryce-carson]

It appears that two separate bugs are in play this round of the great debugging game:

1. The dates of the layers are off-by-one, with the first round of Bayesian data assimilation occurring during the fifth day of August 2018, which is labelled as the sixth day of the month in the pictured plot due to the off-by-one bug. This does not affect when data assimilation actually occurs, the subplots themselves are accurate and the timing is correct, but the post-hoc labelling of the layer's dates is off-by-one.
2. The area of the provinces of Ituri and Nord-Kivu is violated by some procedure in the non-BDA algorithm, causing the (visual) extent to increase. The effect of this on the extent attribute of the SpatRaster object itself is something I will investigate.
3. Bayesian data assimilation appears to be applying the processes of randomization to cells with no value, resulting in them having a value with the pattern of randomization weights and factors exhibited in the Q, H, and HQHT matrices.

With BDA

Psi matrix

lattice::levelplot of the first row of the linear interpolation operator matrix

lattice::levelplot of HQHt

Without BDA

This plot with lines, showing exactly how the values of cells without the extent of the provinces of Ituri and Nord-Kivu are changing due to the second bug relevant to this issue.

[bryce-carson]

The funky pattern of values in the Infected compartment mostly resembles the appearance of the Kalman filter, which has the following appearance.

lattice::levelplot(matrix(KalmanFilter[, 1], byrow = TRUE, ncol = ncol(layers)))

The matrix multiplication of the Kalman filter and the innovation matrix looks like this:

lattice::levelplot(matrix(KalmanFilter %*% Innovation, byrow = TRUE, ncol = ncol(layers)))

It looks more or less identical.

The usage of this is

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 1268 to 1270 in 799d252

	KalmanFilter <- QHt %*% Matrix::solve(HQHt + Psi)
	Posterior <- Prior + KalmanFilter %*% Innovation
	Posterior[Posterior < 0] <- 0

[bryce-carson]

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 1235 to 1245 in 799d252

	assimilateData <-
	function(layers,
	linearInterpolationMatrix,
	prevalenceData,
	healthZoneCoordinates,
	psi.diagonal,
	QHt,
	HQHt) {
	Infected <- terra::as.matrix(layers$Infected, wide = TRUE)
	Prior <- matrix(Matrix::t(Infected), ncol = 1)

The Prior is the row vector of the transposed matrix of values in the Infected layer, which before being linearized looks like this (with the non-BDA bug still affecting the layer).

[bryce-carson]

The issue, then, is located in these lines. I suspect, primarily, line 1269.

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 1246 to 1270 in 799d252

	## FIXME DONE: all elements in the Forecast matrix are NaN. TODO:
	## investigate why the Prior contains NaNs.
	if (any(is.nan(Prior))) {
	warning("Prior contains NaNs, replacing with zeroes.")
	Prior[is.nan(Prior)] <- 0
	}
	Forecast <- linearInterpolationMatrix %*% Prior
	## Create the measurement error covariance matrix.
	Innovation <- as.numeric(prevalenceData) - Forecast
	Psi <- as.numeric(prevalenceData)
	Psi[Psi == 0] <- psi.diagonal
	Psi <- diag(Psi)
	## NOTE: I did not have great documentation left for myself regarding this.
	## Why was I concerned about the dimensions? Why was I concerend whether it
	## was a direct sum or Kronecker sum (or not)?
	if (all(!is.null(getOption("spatialEpisim.foundation.debugMessages")),
	getOption("spatialEpisim.foundation.debugMessages"),
	!identical(dim(HQHt), dim(Psi))))
	message("HQHT + diag(Psi) is either a direct sum or a Kronecker sum, because the dimensions are inequal.")
	KalmanFilter <- QHt %*% Matrix::solve(HQHt + Psi)
	Posterior <- Prior + KalmanFilter %*% Innovation
	Posterior[Posterior < 0] <- 0

[bryce-carson]

Bugged

It appears that two [three] separate bugs are in play this round of the great debugging game:
2. The area of the provinces of Ituri and Nord-Kivu is violated by some procedure in the non-BDA algorithm, causing the (visual) extent to increase. The effect of this on the extent attribute of the SpatRaster object itself is something I will investigate.

Without BDA

This plot with lines, showing exactly how the values of cells without the extent of the provinces of Ituri and Nord-Kivu are changing due to the second bug relevant to this issue.

Fixed (needs to be committed and pushed)

[bryce-carson]

2. The area of the provinces of Ituri and Nord-Kivu is violated by some procedure in the non-BDA algorithm, causing the (visual) extent to increase. The effect of this on the extent attribute of the SpatRaster object itself is something I will investigate.

The cause of this issue was the arguments passed to terra::focal; I did not limit its calculation of values where there were NAs previously. See the function documentation and compare the former and latter code to understand the differences.

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 1011 to 1016 in 0df1abf

	## NOTE: Calculate a matrix of weights respecting human mobility patterns.
	transmissionLikelihoods <-
	terra::focal(x = layers$Infected,
	w = averageEuclideanDistance(lambda, aggregationFactor),
	fun = "sum",
	na.rm = TRUE)

Now it reads,

      ## NOTE: Calculate a matrix of weights respecting human mobility patterns.
      transmissionLikelihoods <-
        terra::focal(x = layers$Infected,
                     w = averageEuclideanDistance(lambda, aggregationFactor),
                     na.policy = "omit",
                     fillvalue = 0,
                     fun = "sum",
                     na.rm = TRUE)