# ---------------------------------------------------------------------
# Program: SurrogatePlot1.S  
#  Author: Steven M. Boker
#    Date: Thu Apr 22 10:32:16 EST 1999
#
# This program creates mutual information plots for the Lorenz
# system and 20 time shuffled surrogates.
# 
#
# ---------------------------------------------------------------------
# Revision History
#   Steven M. Boker -- Thu Apr 22 10:32:12 EST 1999
#      Create SurrogatePlot1.S
#
# ---------------------------------------------------------------------


# -----------------------------
# Read in the The Lorenz attractor

lorenz <- scan("lor63.dat", list(x=0), skip=0)
maxSamples <- length(lorenz$x)
maxSamples
summary(data.frame(x=lorenz$x))

# -----------------------------
# Calculate the Mutual Information (Nonlinear Dependency) for the Lorenz data

maxTau <- 100
write(t(format(round(lorenz$x,6))), file="temp.dat", ncol=1, append=F)
theNLDCommand <- paste("nldcalc -o temp.nld -taumax ", maxTau, " temp.dat", sep="")
unix(theNLDCommand) 
nldLorenz <- scan("temp.nld", list(nld=0), skip=0)
summary(nldLorenz$nld)

# -----------------------------------------------------------------------
# Create 20 time-shuffled surrogates of the Lorenz data
# and calculate their Mutual Information (Nonlinear Dependency)

maxSurrogates <- 20
nldSurrogateMatrix <- matrix(NA, maxSurrogates, (2*maxTau) + 1)

for (i in 1:maxSurrogates) {
	permutedLorenz <- lorenz$x[order(runif(maxSamples))]
	write(t(format(round(permutedLorenz,6))), file="temp.dat", ncol=1, append=F)
	theNLDCommand <- paste("nldcalc -o temp.nld -taumax ", maxTau, " temp.dat", sep="")
	unix(theNLDCommand) 
	nldSurrogateMatrix[i,] <- scan("temp.nld", list(nld=0), skip=0)$nld
}

# -----------------------------
# Plot the Mutual Information (Nonlinear Dependency) for the Lorenz data 
# and time-suffled surrogates

thePlotData <- nldLorenz$nld
tepsfile <- paste("NLDLorenzSurrogates1.eps",sep="")
tmain <- paste("Mutual Information for Lorenz System", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(-maxTau, maxTau), c(0, 0.4), xlab="Tau", ylab="Mutual Information", 
    main=tmain, type="n")
lines(seq(-maxTau, maxTau), thePlotData, type="l")
for (i in 1:maxSurrogates) {
	lines(seq(-maxTau, maxTau), nldSurrogateMatrix[i,], type="l", lty=0)
}
dev.off()

# -----------------------------
# Plot the Mutual Information (Nonlinear Dependency) for the Lorenz data
# and time-suffled surrogates

thePlotData <- nldLorenz$nld
tepsfile <- paste("NLDLorenzSurrogates2.eps",sep="")
tmain <- paste("Mutual Information for Lorenz System", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(-maxTau, maxTau), c(0, 0.05), xlab="Tau", ylab="Mutual Information", 
    main=tmain, type="n")
lines(seq(-maxTau, maxTau), thePlotData, type="l")
for (i in 1:maxSurrogates) {
	lines(seq(-maxTau, maxTau), nldSurrogateMatrix[i,], type="l", lty=0)
}
dev.off()

# -----------------------------
# Calculate the autocorrelation of the Lorenz System

lorenzXX <- rep(NA, maxTau)
for (tau in 1:maxTau) {
	lorenzXX[tau] <- cor(lorenz$x[1:(maxSamples - (tau-1))],
	                     lorenz$x[tau:maxSamples])
}

# -----------------------------
# Paste together the leads and lags

revSeq <- seq(maxTau, 2, -1)
lorenzXX2 <- c(lorenzXX[revSeq],lorenzXX)

# ---------------------------------------------------------------------
# Plot the autocorrelation of the Lorenz system

thePlotData <- lorenzXX2
tepsfile <- paste("LorenzAutoCorr1.eps",sep="")
tmain <- paste("Autocorrelation for Lorenz System", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(-maxTau,maxTau), c(-.1,1), xlab="Tau", ylab="Correlation", main=tmain, type="n")
lines(seq(-(maxTau-1),(maxTau-1)),thePlotData, type="l")
dev.off()


# -----------------------------
# Calculate the autocorrelation of the 20 time shuffled surrogates of 
# the Lorenz System

lorenzXXSurr <- matrix(NA, maxSurrogates, maxTau)
for (i in 1:maxSurrogates) {
	permutedLorenz <- lorenz$x[order(runif(maxSamples))]
	for (tau in 1:maxTau) {
		lorenzXXSurr[i,tau] <- cor(permutedLorenz[1:(maxSamples - (tau-1))],
								permutedLorenz[tau:maxSamples])
	}
}

# -----------------------------
# Paste together the leads and lags of the surrogates

revSeq <- seq(maxTau, 2, -1)
lorenzXXSurr2 <- cbind(lorenzXXSurr[,revSeq],lorenzXXSurr)

# ---------------------------------------------------------------------
# Plot the autocorrelation of the Lorenz system and the 20
# time shuffled surrogates

thePlotData <- lorenzXX2
tepsfile <- paste("LorenzSurrogateAutoCorr1.eps",sep="")
tmain <- paste("Autocorrelation for Lorenz System \n and Time Shuffled Surrogates", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(-maxTau,maxTau), c(-.1,1), xlab="Tau", ylab="Correlation", main=tmain, type="n")
lines(seq(-(maxTau-1),(maxTau-1)),thePlotData, type="l")
for (i in 1:maxSurrogates) {
	lines(seq(-(maxTau-1),(maxTau-1)),lorenzXXSurr2[i,], type="l")
}
dev.off()

thePlotData <- lorenzXX2
tepsfile <- paste("LorenzSurrogateAutoCorr2.eps",sep="")
tmain <- paste("Autocorrelation for Lorenz System \n and Time Shuffled Surrogates", sep="")
postscript(tepsfile,height=6.4,horizontal=F)
plot(c(-maxTau,maxTau), c(-.1,.1), xlab="Tau", ylab="Correlation", main=tmain, type="n")
lines(seq(-(maxTau-1),(maxTau-1)),thePlotData, type="l")
for (i in 1:maxSurrogates) {
	lines(seq(-(maxTau-1),(maxTau-1)),lorenzXXSurr2[i,], type="l")
}
dev.off()


# ---------------------------------------------------------------------
# Quit the program here  

q()