Learn R Programming

pedometrics (version 0.4-1)

spsannMSSD: Optimization of spatial samples for spatial interpolation

Description

This funtion optimizes spatial samples for spatial interpolation using spatial simulated annealing. The criterion used in the optimization is the mean squared shortest distance (MSSD).

Usage

spsannMSSD(points, candidates, x.max, x.min, y.max, y.min, iterations = 10000,
  acceptance = list(initial = 0.99, cooling = iterations/10),
  stopping = list(max.count = iterations/10), plotit = TRUE, boundary,
  progress = TRUE, verbose = TRUE)

objMSSD(candidates, points)

Arguments

points
data.frame or matrix with three columns: 1) the identification of each point, 2) the x coordinates of the set of points, and 3) the y coordinates of the set of points. The coordinates must be projected. This is the set of points to be optimized. poi
candidates
data.frame or matrix with the candidate locations for the sample points. See Details for more information.
x.max,x.min,y.max,y.min
The minimum and maximum quantity of random noise to be added to the x and y coordinates. The minimum quantity should be equal to, at least, the minimum distance between two neighboring candidate locations. The units are the same as of the coordinates. See
iterations
Integer value defining the maximum number of iterations that should be used for the optimization. See Details for more information.
acceptance
List with two sub-arguments: initial and cooling. initial is a numeric value between 0 and 1 defining the initial acceptance probability. Defaults to initial = 0.99. cooling is a numeric val
stopping
A list with one sub-argument: max.count. max.count is an integer value defining the maximum allowable number of iterations without improvement of the objective function value. This is also known as the freezing criterion. Default
plotit
Logical value for ploting the optimization results. This includes a) the progress of the objective function values and acceptance probabilities, and b) the original points, the perturbed points and the progress of the maximum perturbation in the x and y c
boundary
SpatialPolygon defining the boundary of the spatial domain. It is mandatory if plotit = TRUE.
progress
Logical value for printing a progress bar. Defaults to progress = TRUE.
verbose
Logical value for printing messages about the progress of the optimization.

Value

  • objMSSD returns a numeric value: the mean squared shortest distance between a set of points and all grid cells.

    spsannMSSD returns a matrix: the optimized sample points with the evolution of the energy state during the optimization as an attribute.

Spatial simulated annealing

Search graph{ The search graph corresponds to the set of effective candidate locations for a point being jittered in a given iteration. The size of the search graph, i.e. the maximum distance that a point can be moved around, is correlated with the concept of temperature. A larger search graph is equivalent to higher temperatures, which potentially result in more movement or agitation of the set of points or particles.

The current implementation of spatial simulated annealing uses a linear cooling schedule depending upon the number of iterations to control the size of the search graph. The equations are as follows:

x.max.b <- x.max.a - k / iterations * (x.max.a - x.min) y.max.b <- y.max.a - k / iterations * (y.max.a - y.min)

where x.max.a and y.max.a are the maximum allowed shift in the x and y coordinates in the current iteration, x.min and y.min are the minimum required shift in the x and y coordinates, and x.max.b and y.max.b are the maximum allowed shift in the x and y coordinates in the next iteration. iterations is the total number of iterations and k is the current iteration. } Acceptance probability{ The acceptance probability is the chance of accepting a new system configuration that is worse than the current system configuration. The concept of acceptance probability is related with that of temperature. A higher acceptance probability is equivalent to higher temperatures, which potentially result in more movement or agitation of the set of points or particles.

Using a low initial acceptance probability turns the spatial simulated annealing into a greedy algorithm. It will converge in a shorter time, but the solution found is likely to be a local optimum instead of the global optimum. Using a high initial acceptance probability (>0.8) usually is the wisest choice.

An exponential cooling schedule depending upon the number of iterations is used in the current implementation of the spatial simulated annealing to control the acceptance probability. The acceptance probability at each iteration is calculates as follows:

actual_prob <- acceptance$initial * exp(-k / acceptance$cooling)

where actual_prob is the acceptance probability at the k-th iteration, acceptance$initial is the initial acceptance probability, and acceptance$cooling is the exponential cooling factor. } Starting system configuration{ Unidimensional criterion such as the number of points per lag distance class are dependent on the starting system configuration by definition. This means that, depending on the parameters passed to the spatial simulated annealing algorithm, many points will likely to stay close to their starting positions. It would be reasonable to use a starting system configuration that is close to the global optimal, but such thing is not feasible.

Increasing the initial acceptance probability does not guarantee the independence from the starting system configuration. The most efficient option in the current implementation of the spatial simulated annealing algorithm is to start using the entire spatial domain as search graph. This is set using the interval of the x and y coodinates to set x.max and y.max (See above).

An alternative is to start jittering (randomly perturbing) several points at a time and use a cooling schedule to exponentially decrease the number of points jittered at each iteration. The current implementation of the spatial simulated annealing does not explore such alternative. The cooling schedule would be as follows:

new.size <- round(c(old.size - 1) * exp(-k / size.factor) + 1)

where old.size and new.size are the number of points jittered in the previous and next iterations, size.factor is the cooling parameter, and k is the number of the current iteration. The larger the difference between the starting system configuration and the global optimum, the larger the number of points that would need to be jittered in the first iterations. This will usually increase the time spent on the first iterations. } Number of iterations{ The number of iterations has a large influence on the performance of the spatial simulated annealing algorithm. The larger the number of possible system configurations, the higher should the number of iterations be.

The number of possible system configurations increases with:

  • a high initial acceptance probability
  • the use of an infinite set of candidate locations
  • the use of a very dense finite set of candidate locations
}

concept

simulated annealing

Details

Distances{ Euclidean distances between points are calculated. This computation requires the coordinates to be projected. The user is responsible for making sure that this requirement is attained. } Mean squared shortest distance{ Calculating the matrix of distances between all sample points and all prediction locations is computationally expensive. As such, the full matrix of distances is calculated only once for the initial system configuration before the first iteration. At each iteration, only the distance between the new sample point and all prediction locations is calculated. This numeric vector is used to replace the column of the matrix of distances which contained the distances between the old jittered sample point and all prediction locations. The mean squared shortest distance of the new system configuration is then calculated using the updated matrix of distances. The whole proceedure is done at the C-level to speed-up the computation. } Utopia and nadir points{ Knowledge of the utopia and nadir points can help in the construction of multi-objective optimization problems.

the MSSD is a bi-dimensional criterion because it explicitly takes into account both y and x coordinates. It aims at the spread of points in the geographic space. This is completely different from the number of points per lag distance class which is an uni-dimensional criterion -- it aims at the spread on points in the variogram space. It is more difficult to calculate the utopia and nadir points of a bi-dimensional criterion.

The utopia ($f^{\circ}_{i}$) point of MSSD is only known to be larger than zero. It could be approximated using the k-means algorithm, which is much faster than spatial simulated annealing, but does not guarantee to return the true utopia point. The nadir ($f^{max}_{i}$) point is obtained when all sample points are clustered in one of the corners of the spatial domain. This cannot be calculated and has to be approximated by simulation or using the knowledge of the diagonal of the spatial domain (the maximum possible distance between two points).

One alternative strategy is to first optimize a set of sample points using the MSSD as criterion and then create geographic strata. In the multi-objective optimization one would then have to define an unidimensional criterion aiming at matching the optimal solution obtained by the minimization of the MSSD. One such uni-dimensional criterion would be the difference between the expected distribution and the observed distribution of sample points per geographic strata. This criterion would aim at having at least one point per geographic strata -- this is similar to optimizing sample points using the number of points per lag distance class.

A second uni-dimensional criterion would be the difference between the expected MSSD and the observed MSSD. This criterion would aim at having the points coinciding with the optimal solution obtained by the minimization of the MSSD. In both cases the utopia point would be exactly zero ($f^{\circ}_{i} = 0$). The nadir point could be easily calculated for the first uni-dimensional criterion, but not for the second. }

References

Brus, D. J.; de Gruijter, J. J.; van Groenigen, J. W. Designing spatial coverage samples using the k-means clustering algorithm. In: P. Lagacherie, A. M.; Voltz, M. (Eds.) Digital soil mapping - an introductory perspective. Elsevier, v. 31, p. 183-192, 2006.

de Gruijter, J. J.; Brus, D.; Bierkens, M.; Knotters, M. Sampling for natural resource monitoring. Berlin: Springer, p. 332, 2006.

Walvoort, D. J. J.; Brus, D. J.; de Gruijter, J. J. An R package for spatial coverage sampling and random sampling from compact geographical strata by k-means. Computers and Geosciences. v. 36, p. 1261-1267, 2010.

See Also

distanceFromPoints, stratify.

Examples

Run this code
require(sp)
data(meuse.grid)
meuse.grid <- as.matrix(meuse.grid[, 1:2])
meuse.grid <- matrix(cbind(c(1:dim(meuse.grid)[1]), meuse.grid), ncol = 3)
points <- sample(c(1:dim(meuse.grid)[1]), 155)
points <- meuse.grid[points, ]
objMSSD(meuse.grid, points)

Run the code above in your browser using DataLab