# This is an example R source file demonstrating the use of the OUSS process as a null hypothesis for the detection of spectral peaks
# Execute using source("example.R") in the console.
# The actual peak detection and statistical evaluation is done using the peacots package, available on CRAN.
#
# Stilianos Louca
# October 30, 2014

library(peacots)

#Specify example time series (as space-separated values in files)
dataFiles 	= c("Example_data/GPDD_9711", 	"Example_data/GPDD_407", 	"Example_data/cyclic", "Example_data/non_cyclic");

 # specify frequency theshold for low-frequency trimming for each time series
minFreqs	= c(0, 0.03, 0, 0);


#Specify output directory for plots and reports
outputDir = "Example_data_analyzed";
dir.create(outputDir, showWarnings=FALSE);

detrending 	= FALSE; 	# whether to detrend the time series prior to testing

old.par <- par(mfrow=c(length(dataFiles), 3))
for(i in 1:length(dataFiles)){
	#Load particular time series from data file
	filename 	= dataFiles[i];
	timeSeries 	= read.csv(filename, header=FALSE, sep=",", comment.char="#");
	signal_raw 	= timeSeries[,2];
	
	# average time step
	time_step = (tail(timeSeries[,1],n=1)-timeSeries[1,1])/(length(timeSeries[,1])-1);
	
	#detrend if requested. Not recommended for the OUSS test.
	if(detrending){
		fit=loess(Y~X,data.frame(X=timeSeries[,1],Y=timeSeries[,2]),span=0.5,degree=1);
		timeSeries[,2] = fit$residuals;
	}
	
	#Find spectral peak and its statistical significance vs the OUSS process and white noise
	peak_OUSS 	= peacots::evaluate.pm(times=timeSeries[,1], signal=timeSeries[,2], minFitFreq=minFreqs[i], minPeakFreq=minFreqs[i]);
	peak_WN 	= peacots::evaluate.pm.wn(times=timeSeries[,1], signal=timeSeries[,2], minFitFreq=minFreqs[i], minPeakFreq=minFreqs[i]);
	
	#Print results
	cat(sprintf("\n%s:\nSpectral peak detected at frequency %.3g, with power %.3g\nGlobal significance: %.5g vs OUSS, %.5g vs WN. Local significance vs OUSS = %.3g\n  Fitted OUSS: power_o = %.3g, lambda = %.3g, powe_e = %.3g, MLL=%.3g\n  Fitted WN power = %.3g, RSS=%.3g\n", filename, peak_OUSS$frequencies[peak_OUSS$peakMode], peak_OUSS$periodogram[peak_OUSS$peakMode], peak_OUSS$P, peak_WN$P, peak_OUSS$Plocal, peak_OUSS$power_o, peak_OUSS$lambda, peak_OUSS$power_e, peak_OUSS$MLL, peak_WN$powerEstimate, peak_WN$RSS))
	
	
	
	##################
	# Time series
	
	#plot time series
	plot(ts(timeSeries[,1]), ts(timeSeries[,2]), xy.label=FALSE, type="l", ylab="Value", xlab="time", main=sprintf("Time series (%s)",filename), cex=0.8, cex.main=0.9)
	
	
	
	##################
	# OUSS periodogram
	
	#plot periodogram OUSS
	plot(ts(peak_OUSS$frequencies), ts(peak_OUSS$periodogram), xy.label=FALSE, type="l", ylab="Power", xlab="frequency", main=sprintf("Periodogram vs OUSS\n(fmax=%.3g, P=%.3g, Plocal=%.5g)",peak_OUSS$frequencies[peak_OUSS$peakMode],peak_OUSS$P, peak_OUSS$Plocal), col="black", cex=0.8, cex.main=0.9);
	
	#plot fitted OUSS periodogram
	lines(peak_OUSS$frequencies[peak_OUSS$minFitMode:length(peak_OUSS$frequencies)], peak_OUSS$fittedPS[peak_OUSS$minFitMode:length(peak_OUSS$fittedPS)], col="red");
	
	#OUSS plot legend
	legend((0.75*peak_OUSS$frequencies[1]+0.25*tail(peak_OUSS$frequencies,1)), (0.85*max(peak_OUSS$periodogram)), c("periodogram", "fitted OUSS"), lty=c(1,1), col=c("black", "red"), bty="n")
	
	
	
	##################
	# WN periodogram
	
	#plot periodogram WN
	plot(ts(peak_WN$frequencies), ts(peak_WN$periodogram), xy.label=FALSE, type="l", ylab="Power", xlab="frequency", main=sprintf("Periodogram vs WN\n(fmax=%.3g, P=%.3g)",peak_WN$frequencies[peak_WN$peakMode],peak_WN$P), col="black", cex=0.8, cex.main=0.9);
	
	#plot fitted WN periodogram
	lines(c(peak_WN$frequencies[peak_WN$minFitMode],tail(peak_WN$frequencies,1)), c(peak_WN$powerEstimate, peak_WN$powerEstimate ), col="blue");
		
	#WN plot legend
	legend((0.75*peak_WN$frequencies[1]+0.25*tail(peak_WN$frequencies,1)), (0.85*max(peak_WN$periodogram)), c("periodogram", "fitted WN"), lty=c(1,1), col=c("black", "blue"), bty="n")



	##################
	# Save analysis
	
	# write data to file
	dataFile = paste(outputDir, paste0(basename(filename),"_data.txt"), sep="/");
	cat(sprintf("# OUSS and WN analysis of times series %s\n# times\tsignal_raw\tsignal_preprocessed\n", filename), file=dataFile, append=FALSE);
	write.table(data.frame(timeSeries[,1], signal_raw, timeSeries[,2]), , file=dataFile, append=TRUE, row.names=FALSE, col.names=FALSE, quote=FALSE, sep="\t");
	cat(sprintf("\n\n# periodogram of time series used for OUSS test\n# frequency\tpower\tfitted_OUSS_power\n"), file=dataFile, append=TRUE);
	write.table(data.frame(peak_OUSS$frequencies, peak_OUSS$periodogram, ps_ouss(peak_OUSS$frequencies, peak_OUSS$power_o, peak_OUSS$lambda, peak_OUSS$power_e, time_step)), , file=dataFile, append=TRUE, row.names=FALSE, col.names=FALSE, quote=FALSE, sep="\t");
	cat(sprintf("\n\n# periodogram of time series used for WN test\n# frequency\tpower\tfitted_WN_power\n"), file=dataFile, append=TRUE);
	write.table(data.frame(peak_WN$frequencies, peak_WN$periodogram, rep(peak_WN$powerEstimate, length(peak_WN$frequencies))), , file=dataFile, append=TRUE, row.names=FALSE, col.names=FALSE, quote=FALSE, sep="\t");

}

# save overall plot
dev.copy2pdf(file=paste(outputDir, "summary_plot.pdf", sep="/"));



par(old.par)