################################################################
#
#  R script for ENVR workshop Oct., 22. 2008, Boulder, CO
#
################################################################
# Illustration for module by Reinhard Furrer
require( fields)

# first illustration
set.seed(14)
n <- 15
A <- array( runif(n^2),c(n,n)) + diag(n)
A[A < .75] <- 0


A <- as.spam(A)
display( A)


AtA <- t(A) %*% A

display( AtA)

################################################################
# formats
i <- c( 2,4,4,5,5)
j <- c( 1,1,2,1,3)

A <- spam(0,5,5)
A[cbind(j,i)] <- c(.1,.2,.5,.3,.6)
A[cbind(i,j)] <- c(6:10)/10
A <- A+diag.spam(1:5)
print(A)

display(A)

# Classical:
c(A)

# Triplet:
triplet(A,TRUE)

# CSR:
A@rowpointers
A@colindices
A@entries

################################################################
# Timing:

timing <- function(expr)
  as.vector(system.time( for (i in 1:N) expr)[1])


set.seed(14)    
N <- 1000       # how many operations
n <- 999        # matrix dimension
cutoff <- 0.9   # what will be set to 0


A <- array( runif(n^2), c(n,n)) + diag(n)
A[A < cutoff] <- 0
S <- as.spam(A)


timing(A + sqrt(A))
timing(S + sqrt(S))

timing(AtA <- t(A) %*% A)
timing(StS <- t(S) %*% S)

timing(A[1,2] <- .5)
timing(A[n,n-1] <- .5)

timing(S[1,2] <- .5)
timing(S[n,n-1] <- .5)


timing(solve(AtA, rep(1,n)))
timing(solve(StS, rep(1,n)))
timing(chol(AtA))
timing(chol(StS))

display(S)
display(StS)


# an even sparser example...
cutoff <- 0.99
A <- array( runif(n^2), c(n,n)) + diag(n)
A[A < cutoff] <- 0
S <- as.spam(A)


timing(AtA <- t(A) %*% A)
timing(StS <- t(S) %*% S)

timing(UA <- chol(AtA))
timing(US <- chol(StS))

summary(US)

display(StS)
display( as.spam(US))

# the perfect sparse object ...
cutoff <- 0.999
A <- array( runif(n^2), c(n,n)) + diag(n)
A[A < cutoff] <- 0
S <- as.spam(A)

timing(AtA <- t(A) %*% A)
timing(StS <- t(S) %*% S)

timing(UA <- chol(AtA))
timing(US <- chol(StS))

summary(US)

display(StS)
display( as.spam(US))


##########################################################################
# Section about spam

spam.version$version

slotNames("spam")
getSlots("spam")

# We construct our simple example matrix
i <- c( 2,4,4,5,5)
j <- c( 1,1,2,1,3)
A <- spam(0,5,5)
A[cbind(j,i)] <- c(.1,.2,.5,.3,.6)
A[cbind(i,j)] <- c(6:10)/10
A <- A+diag.spam(1:5)

A
slotNames(A)
A@entries
A@colindices
A@rowpointers
A@dimension


################################################################
# Construct example covariance matrices

# Covariance matrix:
n <- 100
x <- seq( 0, 1, l=n )
C <- nearest.dist( x, diag=TRUE, upper=NULL)
C@entries <- Wendland( C@entries, dim=2, k=1)
image(C)


# regular grid precision matrix:
nx <- 5
ny <- 10
b1 <- .1
b2 <- .01
x <- expand.grid(1:nx,1:ny)
C <- diag.spam(nx*ny) + (b1-b2) * nearest.dist(x, delta=1, upper=NULL) + 
         +     b2 * nearest.dist(x, delta=sqrt(2), upper=NULL)
image(C)


# irregular grid precision matrix
n <- 5
incidence <- list( i=c( 2,4,4,5,5), j=c( 1,1,2,1,3), rep(.5,5))
C <- spam( incidence, n, n)
C <- diag.spam(n) + t(C) + C


################################################################
# Illustrate different orderings:

b1 <- .1
b2 <- .01
n <- dim(UScounties.storder)[1]


Aspam <- diag.spam(n) + b1 * UScounties.storder + b2 * UScounties.ndorder

Ummd <- chol(Aspam)
Urcm <- chol(Aspam,pivot="RCM")

summary(Ummd)
summary(Urcm)

display(as.spam( Ummd))
display(as.spam( Urcm))


################################################################
# options in spam:

powerboost()

noquote( format( spam.options()))

# will be implemented in 0.15-3
powerboost <- function(flag="on") {
    if (sys.parent() == 0)
        env <- asNamespace("spam")
    else env <- parent.frame()

    current <- spam.options()
    current[c("safemode","cholsymmetrycheck","cholpivotcheck","eps")] <- 
     if (tolower(flag) %in% c("true","on","an","ein")) {
	list(c(FALSE,FALSE,FALSE),FALSE,FALSE,1e-8)			
   } else { list(c(TRUE,TRUE,TRUE),TRUE,TRUE,.Machine$double.eps) 
  }     
    assign(".Spam", current, envir = env)
 invisible( current)
}

################################################################
# fields example:

precipsubset <- (USprecip[,5]==1)
x <- USprecip[ precipsubset,1:2]
Y <- USprecip[ precipsubset,4]

# We need to eliminate points that are very close
# or set lambda to some positive value like 0.001.
subset <- nearest.dist( x, delta=.05)@colindices
x <- x[-subset,]
Y <- Y[-subset]

out <- mKrig( x,Y, m=1, cov.function="wendland.cov", theta=1.5)


out.p <- predict.surface( out, nx=220, ny=120)
surface(out.p, type='I')
US(add=T)
points(x,pch='.')

# On average (approximately) 
out$nonzero.entries/length(out$residuals)
#        points within the taper range.



# with Krig
out <- Krig( x,Y, m=1, cov.function="wendland.cov", theta=1.5,
            lambda=0)

out.q <- predict.surface( out, nx=220, ny=120)

sum( ( out.q$z-out.p$z)^2, na.rm=T)


N <- 1
timing( out <- mKrig( x,Y, m=1, cov.function="wendland.cov",theta=1.5) )
timing( predict.surface( out, nx=220, ny=120) )
timing( out <- Krig( x,Y, m=1, cov.function="wendland.cov",theta=1.5, lambda=0) )
timing( predict.surface( out, nx=220, ny=120) )





sel <- (x[,1] < -105)&(x[,2] > 37)
xr <- x[sel,]
Yr <- Y[sel]

timing(   out1 <- mKrig( xr,Yr, m=1, theta=1.5) )
timing(   out1.p <- predict.surface( out1, nx=220, ny=120) )
surface(out1.p, type='I')
US(add=T)
points(x,pch='.')

timing(   out2 <- mKrig( xr,Yr, m=1, theta=1.5, cov.function="stationary.taper.cov",
                          Taper.args = list(k=0,theta=3)))
timing(   out2.p <- predict.surface( out2, nx=220, ny=120))
surface(out2.p, type='I')
US(add=T)
points(x,pch='.')


diff <- out1.p
diff$z <- out1.p$z-out2.p$z
surface(diff, type='I')
US(add=T)
points(x,pch='.')




######################################################################
# Reinhard Furrer, 2008-10-01
######################################################################