eye <- function(n,sparse=F)
{
if (sparse)
{
sparseMatrix(i=1:n,j=1:n,x=rep(1,n))
}
else
{
diag(rep(1,n))
}
}
is.symmetric <- function(A)
{
if (nrow(A) == ncol(A))
{
if (any(is.na(A))){
A[is.na(A)] <- 0
}
all(A==t(A))
}
else
{
FALSE
}
}
a.n <- as.numeric
is.square <- function(m) dim(m)[1] == dim(m)[2]
center.conf <- function(XL)
{
llply(XL,function(X) t(scale(t(X),scale=FALSE)))
}
Penalty.M <- function(t , periodic=FALSE)
{
tmp = matrix(0,t,t)
tmp[seq(1,t-1),seq(2,t)] = diag(t-1)
L= -1*diag(t) + tmp
Delta = L %*% L
Delta = Delta[-c(t-1,t),,drop=F]
if (periodic)
{
Delta[,1] <- Delta[,1] + 1
Delta[,t] <- Delta[,t] + 1
}
M <- t(Delta) %*% Delta
M
}
penalty.per.curve.list <- function(DL,XL,params)
{
Dist <- laply(DL, function(D) D)
X <- laply(XL, function(X) X, .drop = FALSE)
pen <- function(Dist,x,params) O(Dist, x , params)
penalties <- aaply( X , 3 , function(x) pen(Dist,array(x,c(params$T,params$D,1)),params))
if (length(rownames(DL[[1]]))>0 ){
penalties <- data.frame(id = rownames(DL[[1]]),penalties=penalties)
} else {
penalties <- data.frame(id = 1:dim(X)[3],penalties=penalties)
}
penalties
}
C.L <- function(DL , XL , params)
{
Dist <- laply(DL, function(D) D)
Dist <- array(Dist,c(params$T,params$N,params$N))
X <- laply(XL, function(X) X, .drop = FALSE)
if (params$T >= 3){
L(Dist , X , params) + params$l * O(Dist , X , params)
} else {
L(Dist, X , params)
}
}
L.L <- function(DL, XL, params)
{
N <- params$N
T <- params$T
Dist <- laply(DL, function(D) D)
Dist <- array(Dist,c(T,N,N))
X <- laply(XL, function(X) X, .drop = FALSE)
D.X <- Dist.matrix(X, params) # TxNxN
stopifnot(dim(D.X) == c(T,N,N))# "Ds dimension error")
sum((D.X - Dist)^2)
}
O.L <- function(DL , XL , params)
{
T <- params$T
D <- params$D
N <- params$N
Dist <- laply(DL, function(D) D)
X <- laply(XL, function(X) X, .drop = FALSE)
#stopifnot( dim(M.DN) == c( T*D , T*D ) )# "M.large dimension error")
#stopifnot( dim(X) == c( T , D , N ) )
n <- dim(X)[3] # if called from f, this is =1 instead of N
if (n == 1)
{
M <- params$M.D
stopifnot( dim(M) == c( T*D , T*D ) )# "M.large dimension error")
}
else
{
M <- params$M.DN
stopifnot( dim(M) == c( T*D*N , T*D*N ) )# "M.large dimension error")
}
Y <- array( X , c( T*D*n,1 ) )
#dim of a.m(Y[,i]): T*Dx1
pen <- t(Y) %*% M %*% Y #length N
#stopifnot( length(pen.per.curve) == n )# "error in calculating penalty of one curve" )
pen
}
L <- function(Dist , X , params)
{
N <- params$N
T <- params$T
D.X <- Dist.matrix( X , params) # TxNxN
stopifnot( dim(D.X) == c(T,N,N) )# "Ds dimension error")
sum( (D.X - Dist)^2 )
}
O <- function(Dist , X , params)
{
T <- params$T
D <- params$D
N <- params$N
#stopifnot( dim(M.DN) == c( T*D , T*D ) )# "M.large dimension error")
#stopifnot( dim(X) == c( T , D , N ) )
n <- dim(X)[3] # if called from f, this is =1 instead of N
if (n == 1)
{
M <- params$M.D
stopifnot( dim(M) == c( T*D , T*D ) )# "M.large dimension error")
}
else
{
M <- params$M.DN
stopifnot( dim(M) == c( T*D*N , T*D*N ) )# "M.large dimension error")
}
Y <- array( X , c( T*D*n,1 ) )
#dim of a.m(Y[,i]): T*Dx1
pen <- t(Y) %*% M %*% Y #length N
#stopifnot( length(pen.per.curve) == n )# "error in calculating penalty of one curve" )
pen
}
Dist.matrix <- function(X,params)
{
T <- params$T
D <- params$D
N <- params$N
stopifnot( dim(X) == c(T , D , N) )
array(aaply(X,1,function(Xt) as.matrix(dist(t(matrix(Xt,D,N))))),c(T,N,N))
}
transform.res <- function(res){
Dist <- laply(res$DL, function(D) D)
X <- laply(res$XL, function(X) X, .drop = FALSE)
D.X <- Dist.matrix(X, res$params)
list(D = Dist, X = X, D.X = D.X)
}
x2list <- function(x,params){
T <- params$T
D <- params$D
N <- params$N
X <- array(x,c(T,D,N))
XL <- alply(X, 1, function(x) matrix(x,D,N))
XL
}
XL2num <- function(XL){
a.n(laply(XL,function(x) as.numeric(x)))
}
Cost <- function(x,DL,params){
XL <- x2list(x,params)
C.L(DL, XL, params)
}
check.local.minimum <- function(res){
x <- XL2num(res$XL)
n <- length(x)
g <- grad(Cost, x, DL = res$DL, params = res$params)
if (all.equal(rep(0,n),g)){
cat("The result is a local minimum of the cost function.\n")
} else {
cat("The automatic test for local minimum failed. \n You are now in a browser to check the actual value of the numerical gradient, g.\n")
browser()
}
}
Dist.List <- function(XL,params){
X <- laply(XL, function(X) X, .drop = FALSE)
D.X <- Dist.matrix(X, params)
DL.X <- alply(D.X,1,function(x) x)
DL.X
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.