Running 'RFsimulate.R' [12s]
  Running 'krige.R' [10s]
  Running 'lgcp.R' [178s]
  Running 'lgm.R' [185s]
  Running 'lgmRaster.R' [80s]
  Running 'likfitLgm.R' [8s]
  Running 'matern.R' [172s]
  Running 'maternGmrfPrec.R' [170s]
  Running 'profLlgm.R' [184s]
  Running 'simLgcp.R' [169s]
  Running 'stackRasterList.R' [6s]
Running the tests in 'tests/lgcp.R' failed.
Complete output:
  > 
  > havePackages = c(
  +     'INLA' = requireNamespace('INLA', quietly=TRUE)
  + )
  > 
  > if(requireNamespace("INLA", quietly=TRUE) ) {
  +   INLA::inla.setOption(num.threads=2)
  +   # not all versions of INLA support blas.num.threads
  +   try(INLA::inla.setOption(blas.num.threads=2), silent=TRUE)
  + }
  > 
  > print(havePackages)
   INLA 
  FALSE 
  > 
  > # as in example
  > require('geostatsp')
  Loading required package: geostatsp
  Loading required package: Matrix
  Loading required package: terra
  terra 1.7.71
  > myPoints = vect(cbind(rbeta(100,2,2), rbeta(100,3,4)))
  > 
  > mycov = rast(matrix(rbinom(100, 1, 0.5), 10, 10), extent=ext(0, 1, 0, 1))
  > names(mycov)="x1"
  > 
  > if(all(havePackages)) {
  +   
  + res = lgcp(
  +     formula=~factor(x1),
  +     data=myPoints, 
  +     grid=20, 
  +     covariates=mycov,
  + 		priorCI=list(sd=c(u=0.1, alpha = 0.01), range=c(0.4, 0.41))
  + )
  + 
  + if(length(res$parameters)) {
  + knitr::kable(res$parameters$summary, digits=3)
  + 
  + plot(res$raster[["predict.exp"]])
  + plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5)
  + }
  + 	
  + # intercept only
  + 
  + res = lgcp(
  +     data=myPoints, 
  +     grid=20, 
  +     covariates=mycov,
  + 		formula=~1,
  + 		priorCI=list(sd=c(0.9, 1.1), range=c(0.4, 0.41))
  + )
  + 
  + if(length(res$parameters)) {
  + 
  + knitr::kable(res$parameters$summary, digits=3)
  + 
  + plot(res$raster[["predict.exp"]])
  + plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5)
  + }
  + }
  > 
  > 
  > # some missing values
  > mycov2 = deepcopy(mycov)
  > temp = values(mycov2)
  > temp[1:4] = NA
  > values(mycov2) = temp
  > 
  > if(all(havePackages)) {
  + 
  + res = lgcp(data=myPoints, grid=20, covariates=mycov2,
  + 		formula=~factor(x1),
  + 		priorCI=list(sd=c(0.9, 1.1), range=c(0.4, 0.41))
  + )
  + 
  + if(length(res$parameters)) {
  + 
  + plot(res$raster[["predict.exp"]])
  + plot(myPoints,add=TRUE,col="#0000FF30",cex=0.5)
  + }
  + }
  > 
  > 
  > data('murder')
  > murder = unwrap(murder)
  > data('torontoPop')
  > torontoPdens = unwrap(torontoPdens)
  > torontoIncome = unwrap(torontoIncome)
  > torontoBorder = unwrap(torontoBorder)
  > 
  > myCov = list(
  +     pop=torontoPdens,
  +     inc = log(torontoIncome)
  + )
  > 
  > formula = ~ inc + offset(pop, log=TRUE)
  > 
  > if(all(havePackages)) {
  + 
  + resL=lgcp(formula, data=murder, 
  +     grid=squareRaster(murder, 30),
  +     covariates=myCov,
  +     border=torontoBorder)
  + 
  + resO = lgcp( ~ inc + pop, 
  +     data=murder, 
  +     grid=squareRaster(murder, 30),
  +     covariates=list(inc=myCov$inc, pop=log(myCov$pop)),
  +     border=torontoBorder)
  + 
  + if(length(resL$parameters)) {
  + 
  + 	rbind(resL$param$summary, resO$param$summary)
  + }
  + 
  + 
  + }	
  > 
  > 
  > # check spdfToBrock
  > 
  > if(requireNamespace('diseasemapping', quietly=TRUE)){
  + 	require('diseasemapping')
  + 	
  + 	data('kentucky')
  + 	kentucky = unwrap(kentucky)
  + 	
  + 	popList = list(
  + 			'2002' = kentucky[,c('M.50', 'M.55', 'M.60')],
  + 			'2005' = kentucky[,c('F.50', 'F.55', 'F.60')]
  + 			)
  + 	for(D in names(popList))
  + 		names(popList[[D]]) = paste('expected_', seq(as.numeric(D)-1, len=3), sep='')
  + 	
  + 	popBrick = spdfToBrick(
  + 			x=popList,
  + 			template=squareRaster(kentucky, 10),
  +     	logSumExpected=FALSE
  + 			)
  + 	sum(popList[['2002']]$expected_2001, na.rm=TRUE)		
  + 	sum(values(popBrick[['expected_2001']]), na.rm=TRUE)*prod(res(popBrick))
  + 	
  + 	popBrick2 = spdfToBrick(
  + 			x=popList,
  + 			template=squareRaster(kentucky, 10),
  +     	logSumExpected=TRUE
  + 	)
  + 	
  + 	sum(unlist(lapply(popList, function(qq) apply(values(qq), 2, sum, na.rm=TRUE))))
  + 	sum(exp(values(popBrick2)), na.rm=TRUE)*prod(res(popBrick2))
  + 			
  + }
  > 	
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     5.79    0.29    6.15 
Running the tests in 'tests/lgm.R' failed.
Complete output:
  > library("geostatsp")
  Loading required package: Matrix
  Loading required package: terra
  terra 1.7.71
  > data("swissRain")
  > swissRain = unwrap(swissRain)
  > swissAltitude= unwrap(swissAltitude)
  > swissLandType = unwrap(swissLandType)
  > 
  > bob = function(x) {
  + 	thepar = x$param
  + 	pdf(tempfile("lgm", tmpdir=".", fileext=".pdf"))
  + 	plot(x$predict[["predict"]], main=
  + 					paste(
  + 							paste(names(thepar), thepar, sep="="),
  + 							collapse=", "),cex.main=0.3
  + 	)
  + 	dev.off()
  + }
  > 
  > 
  > 
  > # specify formula name of raster layer
  > swissFit = lgm(data=swissRain[1:60,], formula=rain~ CHE_alt,
  + 		grid=80, covariates=swissAltitude,
  + 		shape=2,  fixShape=TRUE, 
  + 		boxcox=0.5, fixBoxcox=TRUE, 
  + 		aniso=TRUE)	
  > names(swissFit)
  [1] "predict"    "optim"      "betaHat"    "varParam"   "parameters"
  [6] "data"       "model"      "summary"   
  > swissFit$summary[,1:4]
                        estimate       stdErr      ci0.005      ci0.995
  (Intercept)        5.539667532 1.5186552399  1.627870863 9.451464e+00
  CHE_alt            0.001035103 0.0005719224 -0.000438072 2.508277e-03
  sdNugget           0.968107658           NA  0.641881116 1.460134e+00
  sdSpatial          3.151230623           NA  1.493616604 6.648463e+00
  range/1000        46.171859918           NA 17.631646783 1.209099e+02
  shape              2.000000000           NA           NA           NA
  anisoRatio         8.126224277           NA  2.494065801 2.647706e+01
  anisoAngleRadians  0.642065353           NA  0.505367946 7.787628e-01
  anisoAngleDegrees 36.787634917           NA 28.955450418 4.461982e+01
  boxcox             0.500000000           NA           NA           NA
  > bob(swissFit)
  null device 
            1 
  > 
  > 
  > # specify formula using name of list element
  > 
  > swissFitAgain = lgm(
  +     data=swissRain[1:80,], 
  +     formula=rain~ elev+land,
  + 		grid=80, covariates=list(elev=swissAltitude,land=swissLandType),
  + 		shape=1,  fixShape=TRUE, 
  + 		boxcox=0.5, fixBoxcox=TRUE,
  + 		aniso=TRUE, verbose=FALSE)	
  > names(swissFitAgain)
  [1] "predict"    "optim"      "betaHat"    "varParam"   "parameters"
  [6] "data"       "model"      "summary"   
  > swissFitAgain$summary[,c('estimate','stdErr','Estimated')]
                                              estimate       stdErr Estimated
  (Intercept)                             5.8267433146 1.3533565416      TRUE
  elev                                    0.0006942889 0.0006197647      TRUE
  landCropland/natural vegetation mosaic  0.0000000000           NA     FALSE
  landMixed forests                      -0.5529200263 0.3722759763      TRUE
  landCroplands                          -0.7490923590 0.4604695432      TRUE
  landGrasslands                         -0.7465372004 0.5981679226      TRUE
  landUrban and built-up                 -1.1794725318 0.6416584043      TRUE
  landWater bodies                       -1.2376473219 0.8826197492      TRUE
  landDeciduous broadleaf forest          0.3587127247 0.8997706864      TRUE
  landEvergreen needleleaf forest        -1.0077078508 1.3050293552      TRUE
  landDeciduous needleleaf forest        -2.8208590497 1.2235509390      TRUE
  landPermanent Wetlands                 -3.7422186144 2.0667075839      TRUE
  sdNugget                                0.8932900194           NA      TRUE
  sdSpatial                               2.7000563506           NA      TRUE
  range/1000                             53.4871108877           NA      TRUE
  shape                                   1.0000000000           NA     FALSE
  anisoRatio                              6.6745956225           NA      TRUE
  anisoAngleRadians                       0.6253994861           NA      TRUE
  anisoAngleDegrees                      35.8327510635           NA      TRUE
  boxcox                                  0.5000000000           NA     FALSE
  > bob(swissFitAgain)
  null device 
            1 
  > 
  > swissFitAgain = lgm(
  + 	data=swissRain[1:60,], 
  + 	formula="rain",
  + 		grid=80, covariates=swissAltitude,
  + 		shape=1,  fixShape=TRUE, 
  + 		boxcox=0.5, fixBoxcox=TRUE, 
  + 		aniso=TRUE)	
  > swissFitAgain$summary[,c('estimate','stdErr','Estimated')]
                        estimate       stdErr Estimated
  (Intercept)        5.725580581 1.9988235022      TRUE
  CHE_alt            0.001067468 0.0005765973      TRUE
  sdNugget           0.867983062           NA      TRUE
  sdSpatial          3.593905716           NA      TRUE
  range/1000        71.546074504           NA      TRUE
  shape              1.000000000           NA     FALSE
  anisoRatio         7.981526356           NA      TRUE
  anisoAngleRadians  0.634302759           NA      TRUE
  anisoAngleDegrees 36.342871044           NA      TRUE
  boxcox             0.500000000           NA     FALSE
  > bob(swissFitAgain)
  null device 
            1 
  > 
  > 
  >   
  > 
  > swissFitAgain = lgm(data=swissRain, formula="rain",
  + 		grid=80, covariates=list(elev=swissAltitude,land=swissLandType),
  + 		shape=1,  fixShape=TRUE, 
  + 		boxcox=0.5, fixBoxcox=TRUE, 
  + 		aniso=TRUE)	
  > swissFitAgain$summary[,c('estimate','stdErr','Estimated')]
                                              estimate       stdErr Estimated
  (Intercept)                             6.0525478421 1.1726947666      TRUE
  elev                                   -0.0001028764 0.0005095701      TRUE
  landCropland/natural vegetation mosaic  0.0000000000           NA     FALSE
  landMixed forests                      -0.5140076042 0.3557662719      TRUE
  landGrasslands                         -0.7632888385 0.5420796265      TRUE
  landCroplands                          -0.7377049020 0.4465854844      TRUE
  landUrban and built-up                 -1.3199271082 0.6450922388      TRUE
  landEvergreen needleleaf forest        -1.4648339239 0.7622905027      TRUE
  landWater bodies                       -1.3266755496 0.8774555256      TRUE
  landDeciduous needleleaf forest        -2.3802142458 0.9147492306      TRUE
  landDeciduous broadleaf forest          0.1987287296 0.9072897145      TRUE
  landOpen shrublands                    -1.2190667425 1.2732587969      TRUE
  landPermanent Wetlands                 -4.0866237520 1.7919563808      TRUE
  sdNugget                                0.8691160151           NA      TRUE
  sdSpatial                               2.6680347252           NA      TRUE
  range/1000                             49.1183059982           NA      TRUE
  shape                                   1.0000000000           NA     FALSE
  anisoRatio                              6.2702892780           NA      TRUE
  anisoAngleRadians                       0.6401670343           NA      TRUE
  anisoAngleDegrees                      36.6788692503           NA      TRUE
  boxcox                                  0.5000000000           NA     FALSE
  > bob(swissFitAgain)
  null device 
            1 
  > 
  > 
  > # land type, factor covariate
  > swissRes2 =  lgm(rain ~ elev + factor(land), swissRain, 
  + 		grid=30, 
  + 		covariates=list(elev=swissAltitude,land=swissLandType), 
  + 		boxcox=0.5, fixBoxcox=TRUE, 
  + 		aniso=TRUE
  + )
  > swissRes2$summary
                                               estimate       stdErr      ci0.005
  (Intercept)                              6.0525478421 1.1726947666  3.031886298
  elev                                    -0.0001028764 0.0005095701 -0.001415442
  factor(land)Mixed forests               -0.5140076042 0.3557662719 -1.430400793
  factor(land)Grasslands                  -0.7632888385 0.5420796265 -2.159593425
  factor(land)Croplands                   -0.7377049020 0.4465854844 -1.888032879
  factor(land)Urban and built-up          -1.3199271082 0.6450922388 -2.981574600
  factor(land)Evergreen needleleaf forest -1.4648339239 0.7622905027 -3.428364139
  factor(land)Water bodies                -1.3266755496 0.8774555256 -3.586851205
  factor(land)Deciduous needleleaf forest -2.3802142458 0.9147492306 -4.736452119
  factor(land)Deciduous broadleaf forest   0.1987287296 0.9072897145 -2.138294704
  factor(land)Open shrublands             -1.2190667425 1.2732587969 -4.498764063
  factor(land)Permanent Wetlands          -4.0866237520 1.7919563808 -8.702397508
  sdNugget                                 0.8691160151           NA  0.546757020
  sdSpatial                                2.6680347252           NA  1.509910320
  range/1000                              49.1183059982           NA 20.115531653
  shape                                    1.0000000000           NA           NA
  anisoRatio                               6.2702892780           NA  2.563712758
  anisoAngleRadians                        0.6401670343           NA  0.478839868
  anisoAngleDegrees                       36.6788692503           NA 27.435503505
  boxcox                                   0.5000000000           NA           NA
                                                ci0.995      ci0.025
  (Intercept)                               9.073209386  3.754108335
  elev                                      0.001209689 -0.001101615
  factor(land)Mixed forests                 0.402385584 -1.211296684
  factor(land)Grasslands                    0.633015748 -1.825745383
  factor(land)Croplands                     0.412623075 -1.612996367
  factor(land)Urban and built-up            0.341720384 -2.584284663
  factor(land)Evergreen needleleaf forest   0.498696291 -2.958895855
  factor(land)Water bodies                  0.933500106 -3.046456778
  factor(land)Deciduous needleleaf forest  -0.023976372 -4.173089793
  factor(land)Deciduous broadleaf forest    2.535752163 -1.579526434
  factor(land)Open shrublands               2.060630578 -3.714608127
  factor(land)Permanent Wetlands            0.529150004 -7.598793720
  sdNugget                                  1.381532600  0.610829519
  sdSpatial                                 4.714458336  1.730074601
  range/1000                              119.937569920 24.901820267
  shape                                              NA           NA
  anisoRatio                               15.335777185  3.174956847
  anisoAngleRadians                         0.801494201  0.517412225
  anisoAngleDegrees                        45.922234996 29.645536750
  boxcox                                             NA           NA
                                                ci0.975        ci0.05
  (Intercept)                              8.3509873495  4.1236366020
  elev                                     0.0008958626 -0.0009410446
  factor(land)Mixed forests                0.1832814756 -1.0991910469
  factor(land)Grasslands                   0.2991677062 -1.6549304783
  factor(land)Croplands                    0.1375865634 -1.4722726558
  factor(land)Urban and built-up          -0.0555695534 -2.3810094169
  factor(land)Evergreen needleleaf forest  0.0292280072 -2.7186902221
  factor(land)Water bodies                 0.3931056787 -2.7699614534
  factor(land)Deciduous needleleaf forest -0.5873386990 -3.8848428354
  factor(land)Deciduous broadleaf forest   1.9769838936 -1.2936300481
  factor(land)Open shrublands              1.2764746425 -3.3133910927
  factor(land)Permanent Wetlands          -0.5744537838 -7.0341297043
  sdNugget                                 1.2366177869  0.6464630825
  sdSpatial                                4.1145100275  1.8548577902
  range/1000                              96.8848043347 27.7754603220
  shape                                              NA            NA
  anisoRatio                              12.3833266148  3.5420477207
  anisoAngleRadians                        0.7629218439  0.5371479507
  anisoAngleDegrees                       43.7122017511 30.7763105478
  boxcox                                             NA            NA
                                                 ci0.95         ci0.1
  (Intercept)                              7.9814590822  4.5496790281
  elev                                     0.0007352918 -0.0007559167
  factor(land)Mixed forests                0.0711758384 -0.9699404269
  factor(land)Grasslands                   0.1283528012 -1.4579918325
  factor(land)Croplands                   -0.0031371483 -1.3100272287
  factor(land)Urban and built-up          -0.2588447995 -2.1466460767
  factor(land)Evergreen needleleaf forest -0.2109776257 -2.4417485111
  factor(land)Water bodies                 0.1166103542 -2.4511800522
  factor(land)Deciduous needleleaf forest -0.8755856561 -3.5525125543
  factor(land)Deciduous broadleaf forest   1.6910875072 -0.9640098245
  factor(land)Open shrublands              0.8752576077 -2.8508135470
  factor(land)Permanent Wetlands          -1.1391177998 -6.3831082572
  sdNugget                                 1.1684544226  0.6901337903
  sdSpatial                                3.8377116200  2.0099352869
  range/1000                              86.8611341148 31.5025141130
  shape                                              NA            NA
  anisoRatio                              11.0999429508  4.0182595246
  anisoAngleRadians                        0.7431861180  0.5599019787
  anisoAngleDegrees                       42.5814279528 32.0800203181
  boxcox                                             NA            NA
                                                  ci0.9         pval Estimated
  (Intercept)                              7.5554166561 2.453321e-07      TRUE
  elev                                     0.0005501639 8.400038e-01      TRUE
  factor(land)Mixed forests               -0.0580747815 1.485168e-01      TRUE
  factor(land)Grasslands                  -0.0685858445 1.591088e-01      TRUE
  factor(land)Croplands                   -0.1653825754 9.855935e-02      TRUE
  factor(land)Urban and built-up          -0.4932081396 4.074594e-02      TRUE
  factor(land)Evergreen needleleaf forest -0.4879193367 5.465338e-02      TRUE
  factor(land)Water bodies                -0.2021710470 1.305446e-01      TRUE
  factor(land)Deciduous needleleaf forest -1.2079159373 9.267100e-03      TRUE
  factor(land)Deciduous broadleaf forest   1.3614672836 8.266223e-01      TRUE
  factor(land)Open shrublands              0.4126800621 3.383461e-01      TRUE
  factor(land)Permanent Wetlands          -1.7901392468 2.257583e-02      TRUE
  sdNugget                                 1.0945162495           NA      TRUE
  sdSpatial                                3.5416111857           NA      TRUE
  range/1000                              76.5846171985           NA      TRUE
  shape                                              NA           NA     FALSE
  anisoRatio                               9.7844669785           NA      TRUE
  anisoAngleRadians                        0.7204320900           NA      TRUE
  anisoAngleDegrees                       41.2777181825           NA      TRUE
  boxcox                                             NA           NA     FALSE
  > bob(swissRes2)
  null device 
            1 
  > 
  > 
  >   
  > # simulated data (without a CRS)
  > # and all covariates are in 'data' object
  > myModel = c(intercept=0,variance=2^2,nugget=0.5^2, range=4.5,shape=2, 
  + 		cov1=0.2, cov2=-0.5)
  > covariates = rast(
  + 		xmin=0,ymin=0,xmax=10,ymax=10,
  + 		ncols=200,nrows=200,nlyrs=2, crs='+proj=merc')
  > values(covariates) = c(
  + 		rep(seq(0,1,len=nrow(covariates)), ncol(covariates)),
  + 		rep(seq(0,1,len=nrow(covariates)), 
  + 		rep(nrow(covariates), ncol(covariates))))
  > names(covariates) = c("cov1","cov2")
  > 
  > Npoints = 30
  > set.seed(0)
  > myPoints = vect(cbind(runif(Npoints,0,10), 
  +         seq(0,10, len=Npoints)), crs = '+proj=merc')
  > 
  > # don't test using the randomFields package, it's currently broken on some R builds
  > options(useRandomFields = FALSE)
  > 
  > myPoints = RFsimulate(myModel,myPoints)
  > 
  > values(myPoints) = cbind(
  +    values(myPoints), 
  +   as.data.frame(extract(covariates, myPoints))
  + )
  > 
  > myPoints$y= myModel["intercept"] +
  + 		as.matrix(values(myPoints)[,names(covariates)]) %*% 
  + 		myModel[names(covariates)] +
  + 		myPoints$sim+
  + 		rnorm(length(myPoints), 0, sqrt(myModel["nugget"]))
  > 
  > fitLikfit = likfitLgm(y~cov1+cov2, myPoints,  
  + 		param=c(range=1,nugget=0,shape=1)) 
  > 
  > 
  > 
  > 
  > # run lgm without providing covariates
  > fitMLE =  lgm(
  +     formula=y~ cov1+cov2, 
  +     data=myPoints,
  +     grid=10, covariates=list(), 
  + 		shape=1, fixShape=TRUE)
  > 
  > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"])
  [1] 4.712256 4.712273
  > bob(fitMLE)
  null device 
            1 
  > 
  > 
  > # now give covariates as raster brick
  > fitMLE =  lgm( y~ cov1 + cov2, myPoints, grid=10,  
  + 		covariates=covariates,
  + 		shape=1, fixShape=TRUE)
  > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"])
  [1] 4.712256 4.712273
  > bob(fitMLE)
  null device 
            1 
  > # now give covariates as list
  > fitMLE =  lgm(y~ cov1 + cov2, myPoints, grid=10,   
  + 		covariates=list(cov1=covariates[["cov1"]],
  + 				cov2 = covariates[["cov2"]]),
  + 		shape=1, fixShape=TRUE)
  > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"])
  [1] 4.712256 4.712273
  > bob(fitMLE)
  null device 
            1 
  > 
  > # not remove covariates from data
  > myPoints = vect(crds(myPoints),
  + 		atts=values(myPoints)[,"y",drop=FALSE],
  + 		crs = crs(myPoints))
  > 
  > # now give covariates as raster brick
  > fitMLE =  lgm(y~ cov1 + cov2,  myPoints, grid=10,  
  + 		covariates=covariates,
  + 		shape=1, fixShape=TRUE)
  > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"])
  [1] 4.722069 4.712273
  > bob(fitMLE)
  null device 
            1 
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
    29.37    0.85   30.17 
Running the tests in 'tests/matern.R' failed.
Complete output:
  > library("geostatsp")
  Loading required package: Matrix
  Loading required package: terra
  terra 1.7.71
  > 
  > param = c(range=1, shape=1.5,	anisoRatio=2, anisoAngleDegrees=-25)
  > 
  > matern(c(0, 0.001, 100000), param=param)
  [1] 1.000000 0.999994 0.000000
  attr(,"param")
              range             shape          variance            nugget 
          1.0000000         1.5000000         1.0000000         0.0000000 
         anisoRatio anisoAngleRadians anisoAngleDegrees 
          2.0000000        -0.4363323       -25.0000000 
  > 
  > #x=c(0, 0.001, 100000);param=c(param, variance=1)
  > 
  > #resultFull = .C("matern", as.double(x), as.integer(length(x)),
  > #		as.double(param["range"]), as.double(param["shape"]),
  > #		as.double(param["variance"]))
  > 
  > 
  > # example with raster
  > myraster = rast(nrows=40,ncols=60,xmin=-3,xmax=3,ymin=-2,ymax=2)
  > 
  > # plot correlation of each cell with the origin
  > myMatern = matern(myraster, y=c(0,0), param=param)
  > myMatern[1:3,1:3]
          lyr.1
  1 0.002385719
  2 0.003084067
  3 0.003973211
  4 0.002349589
  5 0.003047536
  6 0.003940191
  7 0.002302637
  8 0.002996124
  9 0.003886871
  >  
  > 
  > 
  > bob = function(x) {
  + thepar = attributes(x)$param
  + pdf(tempfile("matern", tmpdir=".", fileext=".pdf"))
  + plot(x, main=
  + 				paste(
  + 				paste(names(thepar), thepar, sep="="),
  + 				collapse=", "),cex.main=0.5
  + )
  + dev.off()
  + }
  > 
  > bob(myMatern)
  null device 
            1 
  > 
  > 
  > bob(matern(myraster, y=c(1,-0.5), 
  + 				param =	c(range=1, shape=1.5,	anisoRatio=2, anisoAngleDegrees=25)			
  + 						)
  + )
  null device 
            1 
  > 
  > bob(matern(myraster,y= c(0,0), 
  + 				param =	c(range=1, shape=25.1,	anisoRatio=2, anisoAngleDegrees=-25)			
  + 		)
  + )
  null device 
            1 
  > 
  > bob(matern(myraster,y= c(0,0), 
  + 				param =	c(range=0, shape=1.5,	anisoRatio=2, anisoAngleDegrees=-25)			
  + 		)
  + )
  null device 
            1 
  > 
  > bob(matern(myraster, y=c(0,0), 
  + 				param =	c(range=100000, shape=1.5,	anisoRatio=2, anisoAngleDegrees=-25)			
  + 		)
  + )
  null device 
            1 
  > 
  > 
  > 				
  > 				
  > 
  > # correlation matrix for all cells with each other
  > myraster = rast(nrows=4,ncols=6,xmin=-3,xmax=3,ymin=-2,ymax=2)
  > myMatern = matern(myraster, param=c(range=0, shape=2))
  > 
  > dim(myMatern)
  [1] 24 24
  > myMatern[1:3,1:3]
  3 x 3 Matrix of class "dsyMatrix"
       [,1] [,2] [,3]
  [1,]    1    0    0
  [2,]    0    1    0
  [3,]    0    0    1
  > 
  > param = c(range=0.2, shape=1.5)
  > set.seed(0)
  > 
  > mypoints = vect(cbind(runif(10), runif(10)), 
  + 		crs = "epsg:2000",
  + 		atts=data.frame(id=1:10))
  > 
  > myDist = forceSymmetric(as.matrix(distance(mypoints)))
  > 
  > myDist3 = myDist[lower.tri(myDist)]
  > 
  > myMatern1 = matern(mypoints, param)
  > myMatern2 = matern(myDist, param)
  > myMatern3 = matrix(NA, length(mypoints),length(mypoints))
  > myMatern3[lower.tri(myMatern3)] = matern(myDist3, param)
  > 
  > 
  > 
  > class(myMatern1)
  [1] "dsyMatrix"
  attr(,"package")
  [1] "Matrix"
  > class(myMatern2)
  [1] "dsyMatrix"
  attr(,"package")
  [1] "Matrix"
  > class(myMatern3)
  [1] "matrix" "array" 
  > 
  > 
  > (myMatern1 - myMatern2)[1:8, 1:4]
  8 x 4 Matrix of class "dgeMatrix"
                [,1]         [,2]         [,3]         [,4]
  [1,]  0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
  [2,]  0.000000e+00 1.355887e-08 1.491364e-07 1.096928e-08
  [3,]  0.000000e+00 1.491364e-07 9.684880e-02 1.787965e-04
  [4,]  0.000000e+00 1.096928e-08 1.787965e-04 3.393112e-07
  [5,] -4.163336e-17 4.066149e-06 2.550393e-05 2.008934e-06
  [6,]  0.000000e+00 9.514759e-11 2.488279e-04 3.295310e-07
  [7,]  0.000000e+00 7.405549e-07 1.106377e-05 4.935520e-07
  [8,]  0.000000e+00 2.307813e-08 1.273513e-06 1.863947e-08
  > sum((myMatern1 - myMatern2)^2)
  [1] 0.01802854
  > (myMatern1 - myMatern3)[1:8, 1:4]
  8 x 4 Matrix of class "dgeMatrix"
                [,1]         [,2] [,3]          [,4]
  [1,]            NA           NA   NA            NA
  [2,]  0.000000e+00           NA   NA            NA
  [3,]  0.000000e+00 1.665335e-16   NA            NA
  [4,]  0.000000e+00 0.000000e+00    0            NA
  [5,] -4.163336e-17 0.000000e+00    0  0.000000e+00
  [6,]  0.000000e+00 0.000000e+00    0 -6.938894e-18
  [7,]  0.000000e+00 0.000000e+00    0  0.000000e+00
  [8,]  0.000000e+00 0.000000e+00    0  0.000000e+00
  > sum((myMatern1 - myMatern3)^2, na.rm=TRUE)
  [1] 2.951492e-32
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     6.31    0.18    6.51 
Running the tests in 'tests/maternGmrfPrec.R' failed.
Complete output:
  > library('geostatsp')
  Loading required package: Matrix
  Loading required package: terra
  terra 1.7.71
  > matrix(NNmat(7, 7)[,25], 7, 7)
       [,1] [,2] [,3] [,4] [,5] [,6] [,7]
  [1,]    0    0    0    6    0    0    0
  [2,]    0    0    5    4    5    0    0
  [3,]    0    5    3    2    3    5    0
  [4,]    6    4    2    1    2    4    6
  [5,]    0    5    3    2    3    5    0
  [6,]    0    0    5    4    5    0    0
  [7,]    0    0    0    6    0    0    0
  > 
  > myr = squareRaster(ext(0,600,0,300), 60)
  > theNN = NNmat(myr)
  > 
  > 
  > params=c(range = 6*xres(myr),
  + 		cellSize=xres(myr),
  + 		shape=2,
  + 		variance=1600)
  > 
  > 
  > # precision matrix without adjusting for edge effects
  > system.time({precMat = maternGmrfPrec(N=theNN, param=params,
  + 					adjustEdges=FALSE)})
     user  system elapsed 
     0.72    0.04    0.77 
  > 
  > 
  > system.time({theNNadj = NNmat(N=myr, nearest=params['shape']+1, adjustEdges=TRUE)})
     user  system elapsed 
     0.09    0.02    0.11 
  > # and with the adjustment
  > system.time({precMatCorr =maternGmrfPrec(N=theNNadj, param=params, 
  + 					adjustEdges=TRUE)}) 
     user  system elapsed 
     0.99    0.05    1.03 
  > 
  > # find better parameters using a numerical optimizer
  > system.time({precMatAdj =maternGmrfPrec(theNNadj, param=params, 
  + 					adjustEdges='optimalShape') })
     user  system elapsed 
     1.12    0.05    1.17 
  > 
  > attributes(precMat)$info$adjustEdges
  [1] FALSE
  > attributes(precMatAdj)$info$adjustEdges
  [1] "optimalShape"
  > attributes(precMatCorr)$info$adjustEdges
  [1] "theo"
  > 
  > Nx = attributes(theNN)$Nx
  > Ny = attributes(theNN)$Ny
  > midcell = Nx*round(Ny/2) + round(Nx/2) # the middle cell
  > edgecell = Nx*5 + 5 # cell near corner
  > 
  > 
  > # show precision of mid cell
  > precMid=matrix(precMat[,midcell], Ny, Nx, byrow=TRUE)
  > precMid[round(Ny/2)+seq(-3, 3), round(Nx/2)+seq(-3, +3)]
               [,1]          [,2]          [,3]         [,4]          [,5]
  [1,]  0.000000000  0.0000000000  0.0000000000  0.000000000  0.0000000000
  [2,]  0.000000000  0.0000000000  0.0000000000 -0.000125894  0.0000000000
  [3,]  0.000000000  0.0000000000 -0.0003776821  0.001678587 -0.0003776821
  [4,]  0.000000000 -0.0003776821  0.0033571746 -0.008593434  0.0033571746
  [5,] -0.000125894  0.0016785873 -0.0085934344  0.017766776 -0.0085934344
  [6,]  0.000000000 -0.0003776821  0.0033571746 -0.008593434  0.0033571746
  [7,]  0.000000000  0.0000000000 -0.0003776821  0.001678587 -0.0003776821
                [,6]         [,7]
  [1,]  0.0000000000  0.000000000
  [2,]  0.0000000000  0.000000000
  [3,]  0.0000000000  0.000000000
  [4,] -0.0003776821  0.000000000
  [5,]  0.0016785873 -0.000125894
  [6,] -0.0003776821  0.000000000
  [7,]  0.0000000000  0.000000000
  > 
  > if(Sys.info()['user'] =='patrick') {
  + 	
  + # invert to get variance matrices
  + 	precMatAdjInv = solve(precMatAdj)
  + 	precMatCorrInv = solve(precMatCorr)
  + 	precMatInv = solve(precMat)
  + 	
  + # check range of diagonal entries, should be
  + 	params['variance']
  + 	range(diag(precMatInv))
  + 	range(diag(precMatCorrInv))
  + 	range(diag(precMatAdjInv))
  + 	
  + # map marginal variances
  + 	tempA=tempC=tempU = attributes(precMat)$raster
  + 	values(tempC) = diag(precMatCorrInv)
  + 	values(tempU) = diag(precMatInv)
  + 	values(tempA) = diag(precMatAdjInv)
  + 	
  + 	if(!interactive()){
  + 		pdf("maternGmrfMarginalVarRasters.pdf",height=10,width=5)
  + 	}
  + 	par(mfrow=c(3,1))
  + 	plot(tempU)
  + 	plot(tempC)
  + 	plot(tempA)
  + 	if(!interactive()){
  + 		dev.off()
  + 	}
  + 	
  + 	
  + 	
  + # variance matrices
  + 	Ncell = Nx*Ny
  + 	midVec = sparseMatrix(midcell,1,x=1,dims=c(Ncell,1))
  + 	edgeVec = sparseMatrix(edgecell,1,x=1,dims=c(Ncell,1))
  + 	
  + 	therast = attributes(precMat)$raster
  + 	values(therast)=NA
  + 	varRast = c(therast, therast,therast,therast,therast,therast)
  + 	names(varRast) = c('mid','edge','midCor','edgeCor','midAdj','edgeAdj')
  + 	values(varRast[['mid']]) = as.vector(Matrix::solve(precMat, midVec))
  + 	values(varRast[['edge']]) = as.vector(Matrix::solve(precMat, edgeVec))
  + 	values(varRast[['midCor']]) =  as.vector(Matrix::solve(precMatCorr, midVec))
  + 	values(varRast[['edgeCor']]) = as.vector(Matrix::solve(precMatCorr, edgeVec))
  + 	values(varRast[['midAdj']]) =  as.vector(Matrix::solve(precMatAdj, midVec))
  + 	values(varRast[['edgeAdj']]) = as.vector(Matrix::solve(precMatAdj, edgeVec))
  + 	
  + 	
  + 	if(!interactive()){
  + 		pdf("maternGmrfPredRasters.pdf",height=11,width=4)
  + 	}
  + 	par(mfrow=c(7,2))
  + 	
  + 	mMid = matern(varRast[[1]],y=xyFromCell(varRast,midcell),param=params)
  + 	mEdge = matern(varRast[[1]],y=xyFromCell(varRast,edgecell),param=params)
  + 	plot(mMid)
  + 	plot(mEdge)
  + 	
  + 	
  + 	plot(varRast[['mid']])
  + 	plot(varRast[['edge']])
  + 	
  + 	plot(varRast[['mid']]-mMid)
  + 	plot(varRast[['edge']]-mEdge)
  + 	
  + 	plot(varRast[['midCor']])
  + 	plot(varRast[['edgeCor']])
  + 	plot(varRast[['midCor']]-mMid)
  + 	plot(varRast[['edgeCor']]-mEdge)
  + 	
  + 	plot(varRast[['midAdj']])
  + 	plot(varRast[['edgeAdj']])
  + 	plot(varRast[['midAdj']]-mMid)
  + 	plot(varRast[['edgeAdj']]-mEdge)
  + 	
  + 	
  + 	if(!interactive()){
  + 		dev.off()
  + 	}
  + 	
  + 	
  + # compare covariance matrix to the matern
  + 	
  + 	xseqFull = seq(-1.5, 1.5, len=401)*params["range"]
  + 	
  + 	diagdist = sqrt(2)
  + 	
  + 	xymid = xyFromCell(varRast, midcell)
  + 	
  + 	xseq = seq(xmin(varRast)+xres(varRast)/2,xmax(varRast), by=xres(varRast))
  + 	yseq = seq(ymin(varRast)+yres(varRast)/2,ymax(varRast), by=yres(varRast))
  + 	
  + 	vx= extract(varRast, 
  + 			vect(cbind(xseq,
  + 							rep(xymid[2],length(xseq)))
  + 			))
  + 	vy = extract(varRast, 
  + 			vect(cbind(
  + 							rep(xymid[1],length(yseq)),
  + 							yseq)
  + 			)
  + 	)
  + 	
  + 	aseq = seq(0,100,by=xres(varRast))
  + 	aseq = sort(unique(c(aseq, -aseq)))
  + 	
  + 	lur = vect(cbind(xymid[1]+aseq,xymid[2]+aseq))
  + 	lll = vect(cbind(xymid[1]-aseq,xymid[2]+aseq))
  + 	
  + 	vur = extract(varRast, 
  + 			lur		)
  + 	
  + 	
  + 	vll= extract(varRast, 
  + 			lll		)
  + 	xyedg = xyFromCell(varRast, edgecell)
  + 	
  + 	xseqe = seq(xmin(varRast)+xres(varRast)/2,xmax(varRast), by=xres(varRast))
  + 	yseqe = seq(ymin(varRast)+yres(varRast)/2,ymax(varRast), by=yres(varRast))
  + 	
  + 	vxe= extract(varRast, 
  + 			vect(cbind(xseqe,
  + 							rep(xyedg[2],length(xseqe)))
  + 			))
  + 	vye = extract(varRast, 
  + 			vect(cbind(
  + 							rep(xyedg[1],length(yseqe)),
  + 							yseqe)
  + 			)
  + 	)
  + 	
  + 	aseqe = seq(0,100,by=xres(varRast))
  + 	aseqe = sort(unique(c(aseqe, -aseqe)))
  + 	
  + 	lure = vect(cbind(xyedg[1]+aseqe,xyedg[2]+aseqe))
  + 	llle = vect(cbind(xyedg[1]-aseqe,xyedg[2]+aseqe))
  + 	
  + 	vure = extract(varRast, 
  + 			lure		)
  + 	
  + 	vlle= extract(varRast, 
  + 			llle		)
  + 	
  + 	
  + 	
  + 	
  + 	theedg = c('edge','edgeAdj','edgeCor')
  + 	
  + 	
  + 	
  + 	themid = c('mid','midAdj','midCor')
  + 	
  + 	
  + 	thecol =c(mid='red',edge='red',midCor='green',edgeCor='green',
  + 			midAdj='blue',edgeAdj='blue')
  + 	thelwd = c(mid=3,midCor=1,midAdj=1,edge=3,edgeCor=1,edgeAdj=1)
  + 	
  + 	if(!interactive()){
  + 		pdf("maternGmrfPred.pdf",height=5,width=10)
  + 	}
  + 	par(mfrow=c(2,2),mar=c(3,2,0,0))
  + # middle cell, full plot
  + 	plot(xseqFull, matern(xseqFull, param=params),
  + 			type = 'l',ylab='cov', xlab='dist',
  + 			main="matern v gmrf",
  + 			ylim=c(0,params['variance']*1.1),col='grey',lwd=9)
  + 	
  + 	matlines(xseq-xymid[1],vx[,themid],col=thecol[themid],lty=1,
  + 			type='o',lwd=thelwd[themid])
  + 	matlines(yseq-xymid[2],vy[,themid],col=thecol[themid],lty=1,
  + 			type='o',lwd=thelwd[themid])
  + 	
  + 	
  + 	matlines(aseq*diagdist, 
  + 			vur[,themid],col=thecol[themid],lty=2,type='o',
  + 			lwd=thelwd[themid])
  + 	matlines(aseq*diagdist, 
  + 			vll[,themid],col=thecol[themid],lty=2,type='o',lwd=thelwd[themid])
  + 	
  + 	
  + 	
  + 	
  + 	
  + 	
  + 	legend("topright", lty=1, 
  + 			col=thecol[c('edge','edgeCor','edgeAdj')],
  + 			legend=c('vanilla','edge','adj')
  + 	)
  + 	
  + 	
  + # edge cell, full plot
  + 	plot(xseqFull, matern(xseqFull, param=params),
  + 			type = 'l',ylab='cov', xlab='dist',
  + 			main="matern v gmrf",
  + 			ylim=c(0,params['variance']*1.1))
  + 	
  + 	
  + 	
  + 	matlines(xseqe-xyedg[1],vxe[,theedg],col=thecol[theedg],lty=1,
  + 			lwd=thelwd[theedg])
  + 	matlines(yseqe-xyedg[2],vye[,theedg],col=thecol[theedg],lty=1,
  + 			lwd=thelwd[theedg])
  + 	
  + 	
  + 	matlines(aseqe*diagdist, vure[,theedg],col=thecol[theedg],lty=2,
  + 			lwd=thelwd[themid])
  + 	
  + 	matlines(aseqe*diagdist, vlle[,theedg],col=thecol[theedg],lty=2,
  + 			lwd=thelwd[theedg])
  + 	
  + 	# middle cell detail plot
  + 	
  + 	plot(xseqFull, matern(xseqFull, param=params),
  + 			type = 'l',ylab='cov', xlab='dist',
  + 			main="matern v gmrf",
  + 			ylim=params['variance']*c(0.5,1.05),
  + 			xlim=c(-0.5, 0.5)*params['range'],col='grey',lwd=4)
  + 	
  + 	
  + 	
  + 	matlines(xseq-xymid[1],vx[,themid],col=thecol[themid],
  + 			lty=1,
  + 			lwd=thelwd[themid])
  + 	matlines(yseq-xymid[2],vy[,themid],col=thecol[themid],lty=1,
  + 			lwd=thelwd[themid])
  + 	
  + 	
  + 	matlines(aseq*diagdist, vur[,themid],col=thecol[themid],lty=2,
  + 			lwd=thelwd[themid])
  + 	
  + 	matlines(aseq*diagdist, vll[,themid],col=thecol[themid],lty=2,
  + 			lwd=thelwd[themid])
  + 	
  + 	
  + 	
  + 	# edge cells
  + 	
  + 	plot(xseqFull, matern(xseqFull, param=params),
  + 			type = 'l',ylab='cov', xlab='dist',
  + 			main="matern v gmrf",
  + 			ylim=params['variance']*c(0.5,1.05),
  + 			xlim=c(-0.5, 0.5)*params['range'])
  + 	
  + 	
  + 	matlines(xseqe-xyedg[1],vxe[,theedg],col=thecol[theedg],
  + 			lty=1,
  + 			lwd=thelwd[theedg])
  + 	matlines(yseqe-xyedg[2],vye[,theedg],col=thecol[theedg],
  + 			lty=1,
  + 			lwd=thelwd[theedg])
  + 	
  + 	
  + 	matlines(aseqe*diagdist, vure[,theedg],col=thecol[theedg],
  + 			lty=2,
  + 			lwd=thelwd[themid])
  + 	
  + 	matlines(aseqe*diagdist, vlle[,theedg],col=thecol[theedg],
  + 			lty=2,
  + 			lwd=thelwd[theedg])
  + 	
  + 	
  + 	
  + 	if(!interactive()){
  + 		dev.off()	
  + 	}
  + 	
  + # matern variance matrix
  + 	covMatMatern = 
  + 			matern(attributes(precMat)$raster, 
  + 					param=params[c('range','shape','variance')])
  + 	
  + 	covMatMaternAdj = 
  + 			matern(attributes(precMat)$raster, 
  + 					param=attributes(precMatAdj)$info$optimal[
  + 							c('range','shape','variance')])
  + 	
  + 	
  + 	
  + 	prodUncor = covMatMatern %*% (precMat)
  + 	prodCor = covMatMatern %*% (precMatCorr)
  + 	prodAdj = covMatMatern %*% (precMatAdj)
  + 	prodAdjOpt = covMatMaternAdj %*% (precMatAdj)
  + 	
  + 	
  + 	
  + 	thebreaks = c(-100, -0.25, 0.25, 0.75, 1.5, 5, 10, 100)
  + 	
  + 	thecol = terrain.colors(length(thebreaks)-1)
  + 	
  + 	if(!interactive()){
  + 		pdf("maternGmrfDiagRast.pdf",height=8,width=5)
  + 	}
  + 	par(mfrow=c(3,1))
  + 	temp = attributes(precMat)$raster
  + 	values(temp) = diag(prodCor)
  + 	plot(temp,breaks=thebreaks,col=thecol,legend=F)
  + #	mapmisc::legendBreaks('right',breaks=thebreaks, col=thecol)
  + 	
  + 	temp = attributes(precMat)$raster
  + 	values(temp) = diag(prodUncor)
  + 	plot(temp,breaks=thebreaks,col=thecol,legend=FALSE)
  + 	
  + 	temp = attributes(precMat)$raster
  + 	values(temp) = diag(prodAdj)
  + 	plot(temp,breaks=thebreaks,col=thecol,legend=FALSE)
  + 	
  + 	if(!interactive()){
  + 		dev.off()
  + 	}
  + 	
  + 	if(!interactive()){
  + 		pdf("maternGmrfHist.pdf",height=9,width=7)
  + 	}
  + 	par(mfrow=c(4,2),mar=c(4,4,0,0))	
  + 	
  + 	hist(Matrix::diag(prodUncor),breaks=60,ylab='vanilla')
  + 	abline(v=1)
  + 	hist(prodUncor[lower.tri(prodUncor,diag=FALSE)],breaks=60)	
  + 	abline(v=0)
  + 	
  + 	thebreaks = c(-10,seq(0.4, 1.1, len=61),10)
  + 	
  + 	if(interactive()){
  + 		hist(Matrix::diag(prodCor),breaks=thebreaks,ylab='edge')
  + 	} else {
  + 		hist(Matrix::diag(prodCor),breaks=60,ylab='edge')
  + 	}
  + 	abline(v=1)
  + 	hist(prodCor[lower.tri(prodCor,diag=FALSE)],breaks=60)	
  + 	abline(v=0)
  + 	
  + 	if(interactive()){
  + 		hist(Matrix::diag(prodAdj),breaks=thebreaks,ylab='adj')
  + 	} else {
  + 		hist(Matrix::diag(prodAdj),breaks=60,ylab='edge')
  + 	}
  + 	abline(v=1)
  + 	hist(prodAdj[lower.tri(prodAdj,diag=FALSE)],breaks=60)	
  + 	abline(v=0)
  + 	
  + 	if(interactive()){
  + 		hist(Matrix::diag(prodAdjOpt),breaks=thebreaks,ylab='adj optimal')
  + 	} else {
  + 		hist(Matrix::diag(prodAdjOpt),breaks=60,ylab='edge')
  + 	}
  + 	abline(v=1)
  + 	hist(prodAdjOpt[lower.tri(prodAdjOpt,diag=FALSE)],breaks=60)	
  + 	abline(v=0)
  + 	
  + 	
  + 	if(!interactive()){
  + 		dev.off()	
  + 	}
  + 	
  + 	
  + 	# testing marginal variance
  + 	if(FALSE){
  + 		
  + # order one
  + 		Sa = seq(from=4.01,to=5,len=10)
  + 		maxvar=NULL
  + 		for(a in Sa) {
  + 			myp = theNN
  + 			myp@x = c(4+a^2,-2*a,2,1,0,0,0,0)[myp@x]
  + 			myv = solve(myp)
  + 			values(myr) = diag(myv)
  + 			maxvar = c(maxvar,max(myr[15,]))
  + 			
  + 		}
  + 		
  + 		phi = 4/Sa
  + 		scale = sqrt(Sa-4)
  + 		range = sqrt(8)/scale
  + 		
  + 		margvar = (4/pi)^(scale^0.5/2)/(4*pi*(Sa-4))
  + 		margvar = (4/pi)^(scale^2/2)/(4*pi*(Sa-4))
  + 		
  + 		margvarOld = 1/(4*pi*(Sa-4))
  + 		
  + 		
  + 		par(mfrow=c(4,1),mar=c(2,4,0,0))
  + 		they= c(-0.1, 0.1)
  + 		plot(range, maxvar/margvar-1,log='x',ylab='range',ylim=they)
  + 		lines(range, maxvar/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		
  + 		plot(1-phi, maxvar/margvar-1,ylab='phi',ylim=they)
  + 		lines(1-phi, maxvar/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		plot(scale, maxvar/margvar-1,ylab='scale',ylim=they)
  + 		lines(scale, maxvar/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		plot(Sa-4, maxvar/margvar-1,log='x',ylab='Sa-4',ylim=they)
  + 		lines(Sa-4, maxvar/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		max(maxvar/margvar)
  + 		them=which.max(maxvar/margvar)
  + 		c(range[them],phi[them],scale[them],Sa[them])
  + 		
  + 		
  + # order two
  + 		
  + 		Sa2 = seq(from=4.01,to=5,len=40)
  + 		maxvar2=NULL
  + 		for(a in Sa2) {
  + 			myp = theNN
  + 			myp@x = c("1" = a*(a*a+12),
  + 					"2" = -3*(a*a+3),
  + 					"3" = 6*a,
  + 					"4" = 3*a, 
  + 					"5" =  -3,
  + 					"6" = -1)[myp@x]
  + 			myv = solve(myp)
  + 			values(myr) = diag(myv)
  + 			maxvar2 = c(maxvar2,max(myr[15,]))
  + 			
  + 		}
  + 		
  + 		phi2 = 4/Sa2
  + 		scale2 = sqrt(Sa2-4)
  + 		range2 = sqrt(8*2)/scale2
  + 		
  + 		
  + 		margvar = (4/pi)^(scale^2/2)/(4*pi*2*(Sa2-4)^2)
  + 		margvarOld = 1/(4*pi*2*(Sa2-4)^2)
  + 		
  + 		
  + 		par(mfrow=c(4,1),mar=c(2,4,0,0))
  + 		they= c(-0.1, 0.1)
  + 		plot(range, maxvar2/margvar-1,log='x',ylab='range',ylim=they)
  + 		lines(range, maxvar2/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		
  + 		plot(1-phi, maxvar2/margvar-1,ylab='phi',ylim=they)
  + 		lines(1-phi, maxvar2/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		plot(scale, maxvar2/margvar-1,ylab='scale',ylim=they)
  + 		lines(scale, maxvar2/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		plot(Sa-4, maxvar2/margvar-1,log='x',ylab='Sa-4',ylim=they)
  + 		lines(Sa-4, maxvar2/margvarOld-1,col='grey')
  + 		abline(h=0)
  + 		max(maxvar2/margvar)
  + 		them=which.max(maxvar2/margvar)
  + 		c(range2[them],phi2[them],scale2[them],Sa2[them])
  + 		
  + 		
  + 		
  + 		plot(xseqFull, matern(xseqFull, param=params),
  + 				type = 'l',ylab='cov', xlab='dist',
  + 				main="matern v gmrf",
  + 				ylim=params['variance']*c(0.25,1.05),
  + 				xlim=c(-0.75, 0.75)*params['range'])
  + 		
  + 		# middle cell
  + 		themid='mid'
  + 		
  + 		plot(xseq-xymid[1],
  + 				vx[,themid] - matern(xseq-xymid[1], param=params),
  + 				col='red',lty=thelty[themid],
  + 				lwd=thelwd[themid],type='o',pch=1,ylim=c(-55,12))
  + 		matlines(yseq-xymid[2],
  + 				vy[,themid]- matern(yseq-xymid[2], param=params),
  + 				col='orange',lty=thelty[themid],
  + 				lwd=thelwd[themid],type='o',pch=2)
  + 		
  + 		
  + 		matlines(aseq*diagdist, 
  + 				vur[,themid]- matern(aseq*diagdist, param=params),
  + 				col='yellow',lty=thelty[themid],
  + 				lwd=thelwd[themid],type='o',pch=3)
  + 		
  + 		matlines(aseq*diagdist, vll[,themid]- matern(aseq*diagdist, param=params),
  + 				col='pink',lty=thelty[themid],
  + 				lwd=thelwd[themid],type='o',pch=4)
  + 		
  + 		xseqFull = seq(-400,400,len=1000)
  + 		
  + 		mbase =  matern(xseqFull, param=
  + 						c(params[c('variance','range')],
  + 								params['shape'])
  + 		)
  + 		lines(xseqFull, 
  + 				matern(xseqFull, 
  + 						param=
  + 								c(params['variance'],
  + 										params['range']/1.025,
  + 										params['shape']/1.23)
  + 				)-mbase,
  + 				col='blue'
  + 		)
  + 		lines(xseqFull, 
  + 				matern(xseqFull, 
  + 						param=
  + 								c(params['variance'],
  + 										params['range']/1.0225,
  + 										params['shape']/1.22)
  + 				)-mbase,
  + 				col='grey'
  + 		)
  + 	}
  + 	
  + # test providing ar parameter
  + 	
  + 	paramsAR = c(shape=1,oneminusar=0.1)
  + #	N=theNN;param=paramsAR;adjustEdges=TRUE;adjustParam=FALSE 
  + 	
  + 	precMatAR =maternGmrfPrec(theNN, param=paramsAR,
  + 			adjustEdges=FALSE,
  + 			adjustParam=FALSE) 
  + 	
  + 	
  + # find better parameters using a numerical optimizer
  + 	precMatAdjAR =maternGmrfPrec(theNN, param=paramsAR, 
  + 			adjustEdges=TRUE,adjustParam=TRUE) 
  + 	
  + 	
  + # and with the adjustment
  + 	precMatCorrAR =maternGmrfPrec(theNN, param=paramsAR, 
  + 			adjustEdges=TRUE,adjustParam=FALSE) 
  + }
  > 
  > proc.time()
     user  system elapsed 
     9.93    0.43   10.35 
Running the tests in 'tests/profLlgm.R' failed.
Complete output:
  > 
  > library('geostatsp')
  Loading required package: Matrix
  Loading required package: terra
  terra 1.7.71
  > data('swissRain')
  > swissRain = unwrap(swissRain)
  > swissAltitude = unwrap(swissAltitude)
  > 
  > Ncores = c(1,2)[1+(.Platform$OS.type=='unix')]
  > 
  > 
  > 
  > sr2 = swissRain
  > sr2$elev = extract(swissAltitude, sr2)
  Warning message:
  [`[[<-`] only using the first column 
  > swissFit = likfitLgm(
  +     data=sr2, 
  +     formula=rain~ elev,
  +     param=c(range=10000,shape=1,nugget=0,boxcox=0.5,anisoRatio=2,anisoAngleDegrees=45),
  +     paramToEstimate = c("range",'anisoAngleDegrees','anisoRatio'),
  +     reml=FALSE,
  +     verbose=FALSE
  + )
  > 
  > 
  > # calculate log-likelihood at the MLE's, but re-estimate variance
  > sl = loglikLgm(
  +     swissFit$param[c('range','shape','boxcox', 'anisoRatio', 'anisoAngleRadians')],
  +     data=sr2, 
  +     formula=rain~ elev,
  +     reml=swissFit$model$reml)  
  > 
  > 
  > # calculate log-likelihood without re-estimating variance
  > sigSqHat = attributes(sl)$totalVarHat
  > sl1 = loglikLgm(
  +     c(attributes(sl)$param[
  +             c('boxcox','anisoRatio','anisoAngleRadians','shape', 'range')], 
  +         variance=sigSqHat),
  +     data=sr2, 
  +     formula=rain~ elev,
  +     reml=swissFit$model$reml)  
  >  
  > 
  > # re=estimate the anisotropy parameters but not the range
  > sf2 = likfitLgm(
  +     data=swissFit$data, 
  +     formula=swissFit$model$formula,
  +     param= swissFit$param[c('range','nugget','shape','boxcox', 'anisoRatio', 'anisoAngleRadians')],
  +     paramToEstimate = c('variance','anisoAngleRadians','anisoRatio'),
  +     reml=swissFit$model$reml)  
  > 
  > # these should all be the same
  > as.numeric(sl1)
  [1] 644.4812
  > as.numeric(sl) 
  [1] 644.4812
  > swissFit$optim$logL
  m2logL.ml   logL.ml 
   644.4812 -322.2406 
  > sf2$optim$logL
  m2logL.ml   logL.ml 
   644.4812 -322.2406 
  > 
  > date()
  [1] "Sat Mar 23 14:59:45 2024"
  > x=profLlgm(swissFit, mc.cores=Ncores,
  +     range=seq(15000, 55000 , len=12)
  + )
  > date()
  [1] "Sat Mar 23 14:59:49 2024"
  > 
  >  
  > swissInf = informationLgm(swissFit)
  > 
  > 
  > 
  > if(!interactive()) 
  +   pdf("profLswissAngle.pdf")
  > 
  > plot(x[[1]],x[[2]], xlab=names(x)[1],
  + #		yaxt='n',
  + 		ylab='log L',
  + 		ylim=c(min(x[[2]], na.rm=TRUE),x$maxLogL),
  + 		type='n')
  > lines(x[[1]],x[[2]])
  > abline(h=x$breaks[-1],
  + 		col=x$col,
  + 		lwd=1.5)
  > axis(2,at=x$breaks,labels=x$prob,
  + 		line=-1.2,tick=F,
  + 		las=1,padj=1.2,hadj=0,col.axis='red')
  > 
  > abline(v=x$ciLong$par,
  + 		lty=2,
  + 		col=x$col[as.character(x$ciLong$prob)])
  > 
  > 
  > axis(1,at=x$ciLong$par,
  + 		labels=x$ciLong$quantile,
  + 		padj= -6,hadj=0.5, 
  + 		tcl=0.5,cex.axis=0.8,
  + 		col=NA,col.ticks='red',col.axis='red')
  > 
  > ciCols = grep("^ci", colnames(swissInf$summary),
  + 		value=TRUE)
  > 
  > ciValues = unlist(swissInf$summary[intersect(names(x), rownames(swissInf$summary)),ciCols])
  > if(any(!is.na(ciValues)))
  +   axis(1,at=ciValues,
  + 		labels=gsub("^ci","",ciCols),
  + 		padj= 2,hadj=0.5, 
  + 		tcl=-2,cex.axis=0.7,
  + 		col=NA,col.ticks='blue',col.axis='blue')
  > 
  > lines(x[[1]],x[[2]], type='o')
  > 
  > if(!interactive()) 
  +   dev.off()
  null device 
            1 
  > 
  > 
  > if(interactive()  | Sys.info()['user'] =='patrick') {
  +   Ncores = c(1,2)[1+(.Platform$OS.type=='unix')]
  +   date()
  +   x2d=profLlgm(swissFit, mc.cores=Ncores,
  +       anisoAngleRadians=seq(22, 60 , len=24)*(2*pi/360),
  + 		anisoRatio =exp(seq(log(1.5),log(20),len=36))
  +   )
  +   date()
  + if(!interactive()) 
  +   pdf("profLswiss2d.pdf")
  + image(x2d[[1]],x2d[[2]],x2d[[3]],
  + 		breaks=x2d$breaks,
  + 		col=x2d$col,log='y',
  + 		xlab=names(x2d)[1],
  + 		ylab=names(x2d)[2])
  + 
  + thesevars = c("anisoAngleRadians","log(anisoRatio)")
  + thisV = swissInf$information[
  + 		thesevars,thesevars]
  + thisMean= c(x2d$MLE["anisoAngleRadians"],
  + 		log(x2d$MLE['anisoRatio']))
  + 
  + 
  + if(requireNamespace("ellipse", quietly=TRUE)) {
  + 
  + for(D in x2d$prob[x2d$prob>0&x2d$prob<1]) {
  + 	thisE = ellipse::ellipse(thisV, centre=thisMean,
  + 			level=D)
  +   colnames(thisE) = names(thisMean)
  + 	thisE = cbind(thisE,
  + 			anisoRatioExp = exp(thisE[,"anisoRatio"]))
  + 	lines(thisE[,"anisoAngleRadians"],
  + 			thisE[,"anisoRatioExp"],lwd=4)
  + 	lines(thisE[,"anisoAngleRadians"],
  + 			thisE[,"anisoRatioExp"], col=x2d$col[as.character(D)],
  + 			lwd=3)
  + }
  + }
  + 
  + points(x2d$MLE[1],x2d$MLE[2],pch=15) 
  + 
  + legend("topright", fill=x2d$col, legend=x2d$prob[-length(x2d$prob)])
  + 
  + if(!interactive()) 
  +   dev.off()
  + 
  + 
  + 
  + # isotropic
  + sr2$sqrtY = sqrt(sr2$rain)
  + swissFitIso = likfitLgm(
  +     data=sr2, 
  +     formula=sqrtY ~ elev,
  +     param=c(range=10000,shape=1,nugget=0),
  +     reml=FALSE,
  +     verbose=FALSE
  + )
  + 
  + swissFitIso$parameters
  + 
  + newParamList = list(
  + 		range=seq(15, 100 , len=24)*1000,
  +     nugget = seq(0,0.25,len=24)
  + 	)
  + newParam= do.call(expand.grid, newParamList)
  + 
  + res= mapply(
  + 	function(range, nugget) {
  + 		loglikLgm(
  + 		c(swissFitIso$parameters[c('shape')], 
  + 			range = unname(range), nugget = unname(nugget)),
  + 		data = sr2,
  + 		formula = swissFitIso$model$formula,
  + 		reml = swissFitIso$model$reml,
  + 		minustwotimes=FALSE)
  + 	},
  + 	range = newParam$range,
  + 	nugget = newParam$nugget
  + 	)
  + 
  + lMatrix = matrix(res, length(newParamList[[1]]), length(newParamList[[2]]))
  + 
  + if(requireNamespace("mapmisc", quietly=TRUE)) {
  + myCol = mapmisc::colourScale(lMatrix, breaks=8, dec=0)
  + 
  + if(!interactive()) 
  +   pdf("profLswissIso.pdf")
  + image(
  + 	newParamList[[1]]/1000, newParamList[[2]], lMatrix,
  + 	col = myCol$col, breaks=myCol$breaks,
  + 	xlab = names(newParamList)[1],
  + 	ylab = names(newParamList)[2]	
  + 	)
  + mapmisc::legendBreaks('topright', myCol)
  + 
  + if(!interactive()) 
  + 	dev.off()
  + 
  + }
  + }
  > 
  > 
  > proc.time()
     user  system elapsed 
    11.28    0.37   11.68 
Running the tests in 'tests/simLgcp.R' failed.
Complete output:
  > 
  > if(requireNamespace("INLA", quietly=TRUE) ) {
  +   INLA::inla.setOption(num.threads=2)
  +   # not all versions of INLA support blas.num.threads
  +   try(INLA::inla.setOption(blas.num.threads=2), silent=TRUE)
  + } 
  > 
  > library('geostatsp')
  Loading required package: Matrix
  Loading required package: terra
  terra 1.7.71
  > 
  > # exclude this line to use the RandomFields package
  > options(useRandomFields = FALSE)
  > 
  > mymodel = c(mean=-1.5, variance=1, 
  + 				range=2, shape=2)
  > 
  > myraster = rast(nrows=15,ncols=15,xmin=0,xmax=10,ymin=0,ymax=10)
  > 
  > # some covariates, deliberately with a different resolution than myraster
  > covA = covB = myoffset = rast(ext(myraster), 10, 10)
  > values(covA) = as.vector(matrix(1:10, 10, 10))
  > values(covB) = as.vector(matrix(1:10, 10, 10, byrow=TRUE))
  > values(myoffset) = round(seq(-1, 1, len=ncell(myoffset)))
  > 
  > myCovariate = list(a=covA, b=covB, offsetFooBar = myoffset)
  > 
  > set.seed(0)
  > myLgcp=simLgcp(mymodel, myCovariate, 
  +     betas=c(a=-0.1, b=0.25), 
  + 	offset='offsetFooBar',
  + 	rasterTemplate=myraster)
  > 
  > if(requireNamespace("INLA", quietly=TRUE)) {
  + res = lgcp(data=myLgcp$events, 
  + 		formula = ~ a + b + offset(offsetFooBar),
  + 		grid=squareRaster(myoffset, 15), 
  + 		covariates=myCovariate,
  + 		prior=list(sd=0.2, range=0.4))
  + 
  + res$parameters$summary[,c(1,3,5)]
  + 
  + lgcpRoc =  spatialRoc(
  +   fit=res, 
  + 	rr=c(2,3), 
  + 	truth=myLgcp, 
  + 	random=FALSE)
  + 
  + head(lgcpRoc)
  + 
  + plot(lgcpRoc[,'onemspec'] , 
  + 	lgcpRoc[,'2'], 
  + 	type='o', 
  + 	xlim=c(0,1), ylim=c(0,1),
  + 	ylab='sensitivity', xlab='1-specificity'
  + )
  + 
  + }
  > 
  > 
  > 
  > 
  > 
  > 	
  > 
  > proc.time()
     user  system elapsed 
     6.20    0.29    6.46