R Under development (unstable) (2025-08-26 r88710 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> options(error = recover)
> #test = TRUE
> test = FALSE
> mantest = FALSE
> set.seed(1)
> library(intamapInteractive)
Loading required package: intamap
Loading required package: sp
> library(gstat)
> #require(maptools)
> # for SIC2004 dataset
> data(sic2004)
> coordinates(sic.val) = ~x+y
> observations = sic.val["dayx"] 
> coordinates(sic.grid)=~x+y
> predGrid = sic.grid
> 
> #Finding the polygon for the candidate locations
> bb = bbox(predGrid)
> boun = SpatialPoints(data.frame(x=c(bb[1,1],bb[1,2],bb[1,2],bb[1,1],bb[1,1]),
+                                 y=c(bb[2,1],bb[2,1],bb[2,2],bb[2,2],bb[2,1])))
> Srl = Polygons(list(Polygon(boun)),ID = as.character(1))
> candidates = SpatialPolygonsDataFrame(SpatialPolygons(list(Srl)),
+                                       data = data.frame(ID=1))
> 
> # Limits the number of prediction locations to have faster UK 
> # computations
> nGrid = dim(coordinates(predGrid))[1]
> predGrid = predGrid[sample(seq(1,nGrid),1000),]
> # Fits the variogram model (using function fit.variogram from package
> # gstat)
> model = fit.variogram(variogram(dayx~x+y, sic.val), vgm(50, "Sph", 250000, 250))
> #plot(variogram(dayx~x+y, sic.val), model=model)
> # Computes the Mukv of the current network
> initMukv <- calculateMukv(observations, predGrid, model, formulaString = dayx~x+y)
> print(initMukv)
[1] 116.712
> # Computes Kriging and plot result 
> #GammaDoseMap = krige(dayx~x+y, observations, predGrid, model)
> #GammaDoseMap = as.data.frame(GammaDoseMap)
> #windows()
> #levelplot(var1.pred~x+y, GammaDoseMap, aspect = "iso", col.regions=bpy.colors,
> #               panel = function(...) {
> #                       panel.levelplot(...)
> #                       panel.xyplot(y=sic.val$y, x=sic.val$x, col="white", pch=19);  panel.xyplot(y=sic.val$y, x=sic.val$x, col="black")
> #               }, main = "universal kriging prediction")
> #windows()
> #levelplot(var1.var~x+y, GammaDoseMap, aspect = "iso",
> #col.regions=bpy.colors,
> #               panel = function(...) {
> #                       panel.levelplot(...)
> #                       panel.xyplot(y=sic.val$y, x=sic.val$x, col="white", pch=19);  panel.xyplot(y=sic.val$y, x=sic.val$x, col="black")
> #                       
> #               }, main = "universal kriging variance")
> 
> 
> 
> 
> 
> 
> ###############################################################
> # Deleting
> ###############################################################
> 
> # Deletes manually 20 stations from current network with method
> # "manual" 
> if (mantest) {
+ optimDel1=optimizeNetwork( observations,
+                            method = "manual",
+                            action = "del",
+                            nDiff = 2,
+                            predGrid, candidates, plotOptim = FALSE,  formulaString = dayx~x+y)
+ # Computes the Mukv of the optimized network with spatial # coverage
+ MukvDel1 <- calculateMukv(optimDel1, predGrid, model, formulaString = dayx~x+y)
+ print(MukvDel1)
+ }
> 
> # Deletes optimally 20 stations from current network with method
> # "spcov" (spatial coverage)
> #if (TRUE) {
> optim1 = optimizeNetwork(observations,
+          	                method = "spcov",
+                           action = "del",
+                           nDiff = 2,
+                           predGrid, candidates,
+                           plotOptim=FALSE)
> # Computes the Mukv of the optimized network with spatial coverage
> MukvDel2 = calculateMukv(optim1, predGrid, model, formulaString = dayx~x+y)
> print(MukvDel2)
[1] 116.8215
> 
> # Deletes optimally 20 stations from current network with method "ssa"
> # (spatial simulated annealing) and criterion "mukv"
> #windows()
> if (test) {
+ optim2 = optimizeNetwork(observations ,
+                           method = "ssa",
+                           criterion = "MUKV",
+                           action = "del",
+                           nDiff = 2,
+                           predGrid, candidates, model,
+                           plotOptim=FALSE)
+ # Computes the Mukv of the optimized network with spatial simulated
+ # annealing applied to mukv
+ MukvDel3 <- calculateMukv(optim2, predGrid, model)
+ print(MukvDel3)
+ }
> ###############################################################
> # Adding
> ###############################################################
> 
> # Adds manually 20 stations from current network with method
> # "manual" 
> 
> if (mantest) {
+ optimAdd1=optimizeNetwork( observations,
+                            method = "manual",
+                            action = "add",
+                            nDiff = 2,
+                            predGrid, candidates)
+ # Computes the Mukv of the optimized network with spatial # coverage
+ MukvAdd1 <- calculateMukv(optimAdd1, predGrid, model)
+ print(MukvAdd1)
+ }
> # Adds optimally 20 stations from current network with
> # method "spcov" (spatial coverage)
> 
> if (TRUE) {
+   optimAdd2=optimizeNetwork( observations,
+                            method = "spcov",
+                            action = "add",
+                            nDiff = 2,
+                            predGrid, candidates, nGridCells = 5000,
+                            nTry = 100, plotOptim=FALSE)
+ # Computes the Mukv of the optimized network with spatial # coverage
+ MukvAdd2 <- calculateMukv(optimAdd2, predGrid, model)
+ print(MukvAdd2)
+ }
[1] 116.0092
Warning messages:
1: 1 location(s) outside the target universe (as defined by 'object') have been found
These locations have been removed. 199 location(s) have been retained 
2: (Nearly) coinciding points have been removed.  
> 
> # Adds optimally 20 stations from current network with
> # method "ssa" (spatial simulated annealing) and
> # criterion "mukv"
> if (test) {
+ optimAdd3=optimizeNetwork( observations ,
+                            method = "ssa",
+                            criterion = "MUKV",
+                            action = "add",
+                            nDiff = 2,
+                            predGrid, candidates, model,
+                            plotOptim=FALSE, nr_iterations = )
+ 
+ # Computes the Mukv of the optimized network with spatial
+ # simulated annealing applied to mukv
+ MukvAdd3 <- calculateMukv(optimAdd3, predGrid, model)
+ print(MukvAdd3)
+ }
> 
> # Compares computed designs based on the Mukv results (Small MUKV is better)
> print(initMukv)
[1] 116.712
> # Deleting 20 measurements
> #print(MukvDel1) # Manual
> print(MukvDel2) # Spatial Coverage
[1] 116.8215
> if (test) print(MukvDel3) # Spatial simulated annealing (MUKV in objective function)
> # Adding 20 measurements
> #print(MukvAdd1) # Manual
> print(MukvAdd2) # Spatial Coverage
[1] 116.0092
> if (test) print(MukvAdd3) # Spatial simulated annealing (MUKV in objective function)
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   7.35    0.65    7.50