R/libEmdCwa.R

In martindurocher/emdCwa: Hybrid method for analyzing teleconnection

newEmdCwa <- function(tt, index, site)
{
	if(length(index) != nrow(site))
	{ 
		warning('Bad dimension of site dataset')
		return(NULL)
	}
	
	if(length(index) != length(tt)) 
	{ 
		warning('Bad dimension of index')
		return(NULL)
	}
	
	if(sum(is.na(index))!= 0)
	{ 
		warning('NA values in index')
		return(NULL)
	}
	
	if(sum(is.na(site))!= 0)
	{ 
		warning('NA values in site')
		return(NULL)
	}
	
	ans <- list(tt = tt, 
		index = as.numeric(index),
		site = as.matrix(site))
		
	class(ans) <- append(class(ans),'emdcwa')
	
	return(ans)
}

#x : time series
#p : number of imf to extract
#n : number of 
#ratio : ratio of standard deviation in noise
#track : if does print a progress bar
eemd <- function(x,...) UseMethod('eemd',x)
eemd.numeric <- function(x, p, k = 5, ratio = 0.1, track = FALSE, ...)
{
	xSd <- sd(x)
	xLen <- length(x)
	
	ans <- dMat<- matrix(0,xLen,p)
	
	if(track) bar <- txtProgressBar(0,k) 
	
	for(i in seq(k))
	{
		noise <- rnorm(xLen,0, ratio * xSd)
		decomp <- emd(x+noise, max.imf = p-1, ...)

		tmp <- cbind(decomp$imf,decomp$residue)
		dId <- seq(pmin(p,ncol(tmp)))
		dMat[,dId] <- tmp[,dId]
		
		ans <- ans + dMat
		
		if(track) setTxtProgressBar(bar,i)
	} 
	return(ans/k)
}

eemd.emdcwa <- function(self, p, k = 5, ratio = 0.1, track = TRUE, ...)
{
	bar <- txtProgressBar(0,ncol(rain)+1)
	
	naoEmd <- eemd(self$index, p = p, k = k, ...)
	
	setTxtProgressBar(bar,1)
	
	siteEmd <- array(NA,dim = c(ncol(self$site),nrow(naoEmd),ncol(naoEmd)))
	
	for(i in seq(ncol(self$site)))
	{
		siteEmd[i,,]<- eemd(self$site[,i], p = p, k = k, ...)
		setTxtProgressBar(bar,i+1)
	}

	return(setEmd(self, naoEmd,siteEmd))
}

setEmd <- function(self,indexEmd, siteEmd)
{
	self$indexEmd <- as.matrix(indexEmd)
	self$siteEmd <- as.array(siteEmd)
	return(self)
}

cwt <- function(self, m = 1,  kSite = 0, kIndex=0, ...)
{
	if(kSite == 0) 
		ySite <- cbind(self$tt, self$site[,m])
	else
		ySite <- cbind(self$tt, self$siteEmd[m,,kSite])
		
	if(kIndex == 0) 
		yIndex <- cbind(self$tt, self$index)
	else
		yIndex <- cbind(self$tt, self$indexEmd[,kIndex])
	
	return(xwt(ySite,yIndex, ...))
}

multiCwt <- function(self, kSite, kIndex = kSite)
{
	spec <- cwt(self,m=1, kSite = kSite, kIndex = kIndex)
	wtMat <- matrix(NA,ncol(self$site),prod(dim(spec$power.corr)))
	sigma <- matrix(NA,ncol(self$site),2)
	wtMat[1,] <- c(spec$power.corr)
	sigma[1,1] <- spec$d1.sigma
	sigma[1,2] <- spec$d2.sigma
	
	for(i in seq(2,ncol(self$site)))
	{
		spec <- cwt(self,m=i, kSite = kSite, kIndex = kIndex)
		wtMat[i,] <- c(spec$power.corr)
		sigma[i,1] <- spec$d1.sigma
		sigma[i,2] <- spec$d2.sigma	
	}
		
	ans <- list(wtMat = wtMat, tt = self$tt, period = spec$period,
		sigma = sigma)
	
	class(ans) <- append(class(ans),'mcwt')
	return(ans)
}

plotWt <-function(wt, x, y, sigma1, sigma2)
{
	figPer <- c(3,4,6,8,11,16)
	figAt <- log2(figPer)
	figY <- log2(y)
	yLim <- log2(c(max(y),min(y)))
	
	wtMat <- t(matrix(wt,length(y),length(x)))
	wtMat <- wtMat / (sigma1 * sigma2)
	
	#par(mar = c(5,6,3,3.5), oma =c(0,0,0,1))
	image(x = x, y = figY, z = wtMat, ylim = yLim, col = tim.colors(32),
		xlab = 'Time', ylab = "Period (year)", axes = FALSE)
	box()
	axis(side = 1)
	axis(side = 2, at = figAt, labels = figPer)
	image.plot(x= x, z=wtMat, legend.only=TRUE)
	
	print(yLim)
}

prcomp.mcwt <- function(self, ...)
{
	
	return(0)
}

print.pcMcwt <- function(self) return(0)

plot.pcMcwt <- function(self) return(0)



cancor.mcwt <- function(self, ...)	return(0)
print.ccaMcwt <- function(self) return(0)
plot.ccaMcwt <- function(self) return(0)


hclust.mcwt <- function(self) return(0)
print.pcMcwt <- function(self) return(0)
plot.pcMcwt <- function(self) return(0)


overlayPlot <- function(self, m , k, ylim = NULL, ylab = NULL,
	xlab = NULL, ...)
{
	if(is.null(xlab)) xlab <- 'Time' 
	
	if(k == 0)
	{
		if(is.null(ylim))
		{
			vec <- c(self$index,self$site[,m])
			mn <- min(vec)
			mx <- max(vec)
			ylim <- c(mn-0.05*abs(mn), mx+0.05*abs(mx))
		}
				 
		plot(self$tt, self$index, ylim = ylim, xlab = xlab, 
			ylab =  ylab, ...)
		lines(self$tt,self$site[,m])
	
	}
	else
	{
		if(is.null(ylim))
		{
			vec <- c(self$indexEmd[,k],self$siteEmd[m,,k])
			mn <- min(vec)
			mx <- max(vec)
			ylim <- c(mn-0.05*abs(mn), mx+0.05*abs(mx))
		}
		
		plot(self$tt, self$indexEmd[,k],  xlab = xlab, 
			ylab =  ylab, ...)
			
		lines(self$tt,self$siteEmd[m,,k])
	}
	
	return(0)
}

waveletPlot <- function(self, type = 'site', m = 1)
{
	return(0)
}