spsannPPL(points, candidates, lags = 7, lags.type = "exponential",
lags.base = 2, cutoff = NULL, criterion = "distribution",
pre.distri = NULL, 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)objPoints(points, lags = 7, lags.type = "exponential", lags.base = 2,
cutoff = NULL, criterion = "distribution", pre.distri = NULL)
pointsPerLag(points, lags = 7, lags.type = "exponential", lags.base = 2,
cutoff = NULL)
poilags = 7"equidistant", for equidistant lag distance
classes, and "exponential", for exponentially spaced lag distance
classes. Defaults to lags.type =lags.base = 2. See "minimum" and "distribution". The first
returns the minimum number of points or criterion = "distribution". Defaults to a uniforminitial 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 valmax.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. Defaultplotit = TRUE.progress = TRUE.pointsPerLag and pairsPerLag return a data.frame with three
columns: a) the lower and b) upper limits of each lag distance class, and
c) the number of points or point-pairs per lag distance class.objPoints and objPairs return a numeric value depending on the
choice of criterion. If criterion = "distribution", the sum of
the differences between the pre-specified and observed distribution of counts
of points or point-pairs per lag distance class. If
criterion = "minimum", the inverse of the minimum count of points or
point pairs over all lag distance classes multiplied by a constant.
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:
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.
}
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:
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.
}
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:
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.
}
The number of possible system configurations increases with:
points, and $lag$ is the number of lag distance
classes.Using the default uniform distribution means that the number of
points per lag distance class is equal to the total number of points
in points. This is the same as expecting that each point contributes
to every lag distance class.
Distributions other that the available options can be easily implemented
changing the arguments lags, lags.base and pre.distri.
}
lags.type = "equidistant") are evenly spaced lags. They are created
by simply dividing the distance interval from 0.0001 to cutoff by the
required number of lags. The minimum value of 0.0001 guarantees that a point
does not form a pair with itself.
The second type (lags.type = "exponential") of lag distance classes
is defined by exponential spacings. The spacings are defined by the base
$b$ of the exponential expression $b^n$, where $n$ is the
required number of lags. The base is defined using argument lags.base.
For example, the default lags.base = 2 creates lags that are
sequentially defined as half of the immediately preceding larger lag. If
cutoff = 100 and lags = 4, the upper limits of the lag distance
classes will be
objPoints and objPairs (to be implemented) and
were designed to be used in spatial simulated annealing to optimize spatial
sample configurations. Both of them have two criteria implemented. The first
is called using criterion = "distribution" and is used to minimize the
sum of differences between a pre-specified distribution and the observed
distribution of points or point-pairs per lag distance class.
Consider that we aim at having the following distribution of points per lag:
desired <- c(10, 10, 10, 10, 10),
and that the observed distribution of points per lag is the following:
observed <- c(1, 2, 5, 10, 10).
The objective at each iteration of the optimization will be to match the two distributions. This criterion is of the same type as the one proposed by Warrick and Myers (1987).
The second criterion is called using criterion = "minimum". It
corresponds to maximizing the minimum number of points or point-pairs
observed over all lag distance classes. Consider we observe the following
distribution of points per lag in the first iteration:
observed <- c(1, 2, 5, 10, 10).
The objective in the next iteration will be to increase the number of points in the first lag ($n = 1$). Consider we then have the following resulting distribution:
resulting <- c(5, 2, 5, 10, 10).
Now the objective will be to increse the number of points in the second lag ($n = 2$). The optimization continues until it is not possible to increase the number of points in any of the lags, that is, when:
distribution <- c(10, 10, 10, 10, 10).
This shows that the result of using criterion = "minimum" is similar
to using criterion = "distribution". However, the resulting sample
pattern can be significantly different. The running time of each iteration
can be a bit longer when using criterion = "distribution", but since
it is a more sensitive criteria (it takes all lags into account), convergence
is likely to be attained with a smaller number of iterations. Note that this
also depends on the other parameters passed to the optimization algorithm.
It is important to note that using the first criterion
("distribution") in simulated annealing corresponds to a
minimization problem. On the other hand, using the second criterion
("minimum") would correspond to a maximization problem. We
solve this inconsistency substituting the criterion that has to be maximized
by its inverse. For conveninence we multiply the resulting value by a
constant (i.e. $c / x + 1$, where c is the number of points and
x is the criterion value). This procedure allows us to define both
problems as minimization problems.
}
When criterion = "distribution", the utopia
($f^{\circ}_{i}$) point is exactly zero ($f^{\circ}_{i} = 0$). When
criterion = "minimum", the utopia point is approximately 1 (0.9)
($f^{\circ}_{i} \sim 1$). It can be calculated using the equation
$n / n + 1$, where n is the number of points
(objPoints), or the number point-pairs divided by the number of lag
distance classes (objPairs).
THIS SECTION ABOUT THE NADIR POINT HAS TO BE REVIEWED!
The nadir ($f^{max}_{i}$) point depends on a series of
elements. For instance, when criterion = "distribution", if the
desired distribution of point or point-pairs per lag distance class is
pre.distribution <- c(10, 10, 10, 10, 10), the worst case scenario
would be to have all points or point-pairs in a single lag distance class,
that is, obs.distribution <- c(0, 50, 0, 0, 0). In this case, the
nadir point is equal to the sum of the differences between the two
distributions:
sum((c(10, 10, 10, 10, 10) - c(0, 50, 0, 0, 0)) ^ 2) = 2000.
When objective = "minimum", the nadir point is equal to
$f^{max}_{i} = n / 0 + 1 = n$.
}
Marler, R. T.; Arora, J. S. Function-transformation methods for multi-objective optimization. Engineering Optimization. v. 37, p. 551-570, 2005.
Russo, D. Design of an optimal sampling network for estimating the variogram. Soil Science Society of America Journal. v. 48, p. 708-716, 1984.
Truong, P. N.; Heuvelink, G. B. M.; Gosling, J. P. Web-based tool for expert elicitation of the variogram. Computers and Geosciences. v. 51, p. 390-399, 2013.
Warrick, A. W.; Myers, D. E. Optimization of sampling locations for variogram calculations. Water Resources Research. v. 23, p. 496-500, 1987.
dist.require(sp)
data(meuse)
meuse <- as.matrix(meuse[, 1:2])
meuse <- matrix(cbind(c(1:dim(meuse)[1]), meuse), ncol = 3)
pointsPerLag(meuse, cutoff = 1000)
objPoints(meuse, cutoff = 1000)Run the code above in your browser using DataLab