#Written by Vadim Teverovsky

# Contains the top level functions plotfract and plotfract2.
# Both implement Higuchi's algorithm for finding H using the fractal
# dimension.   plotfract (and Sfract) take a non-stationary series of the
# form, e.g., of  Fractional Brownian Motion, and give the H corresponding
# to its difference, e.g. Fractional Gaussian Noise.  plotfract2 (and
# Sfract2) take a stationary series like FGN, or FARIMA.  For long series, if
# the FBM is available, plotfract is much faster, but if only the stationary
# series is available, it is easier to use plotfract2.

plotfract_function(data, nvar = 50, power1 = 0.7, power2 = 2.5, 
		slopes = 1, plotflag = T)
{

#	This is a function which gives a plot of Higuchi's fractal
#	dimension algorithm for various
#	aggregation levels.  It calls the function Sfract.S, which, in
#	turn, calls Cfract.c.  The estimate of the slope, beta, gives
#	an estimate for the parameter, H, via :  H = 2 + beta
#
#NOTE::::   Data should not be of FGN type, but FBM(Fractional Brownian Motion)
#
#	Variables:
#		data	--	The data being used.
#		nvar	-- 	The number of aggregation levels.
#		power1,power2-- Cut-offs for points to be used in estimating
#				beta.
#		slopes  --	Flag.  If 1, draw the estimated line whose 
#				slope is beta.
#		plotflag   --	Flag.  If T, plots the estimates vs. aggregation
#				level.
	cat("\n")	

# n is the length of the data vector, absolute is the vector of length nvar+1
# which is the output of Sfract.  It is the vector of curve lengths  at the 
# different levels of aggregation.
	n <- length(data)
	absolute <- Sfract(data, nvar)	

# The increment to be used in plotting and estimating the slope.  Depends on 
# the number of levels.
	increment <- (log10(n/2))/nvar	

# lv and lm are going to be the y and x values respectively of the points in the
# plot.
	lv <- log10(absolute)
	lm <- rep(0, nvar)	

# power1 and power2 determine whether the points will be included in the 
# estimate of beta.
	for(i in 1:nvar)
		lm[i] <- log10(floor(10^(i * increment)))

#Above line changed 2/28/95, to make the plotting consistent with calculations.

	lmadj <- lm[(lm >= power1) & (lm <= power2)]
	lvadj <- lv[(lm >= power1) & (lm <= power2)]
	lmadjc <- lm[(lm < power1) | (lm > power2)]
	lvadjc <- lv[(lm < power1) | (lm > power2)]	

# The line is fitted to the appropriate points.
	fit <- lsfit(lmadj, lvadj)
	intercept <- fit$coef[1]
	beta <- fit$coef[2]	

#Plot only if plotflag == T.  VT. 12/28/94.

	if(plotflag == T) {
		plot(c(0, lm), c(0, lv - lv[1]), xlab = "log10(m)", ylab = 
			"log10(Curve Length)", type = "n")
		points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.06)
		points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06)
		abline(0, -1.5, lty = 2)	

#The least squares fitted line is drawn if slopes = 1.

		if(slopes == 1)
			lines(lmadj, intercept - lv[1] + (beta) * lmadj)
	}

#Output results to screen.
	cat(" beta = ", beta, "\n\n")
	cat("H = ", 2 + beta, "\n\n")	

#returning a value for the function.VT. 12/28/94
	result <- beta + 2
	names(result) <- NULL
	return(result)
}


Sfract <- function(data, nvar)
{
	n <- length(data)
	if(!is.loaded(C.symbol("Cfract")))
	dyn.load("Cfract.o")
	.C("Cfractal",
		data = as.double(data),
		n = n,
		nvar = as.integer(nvar),
		v = double(nvar + 1))$v[2:(nvar+1)]
}




###########################

plotfract2_function(data, nvar = 50, power1 = 0.7, power2 = 2.5, slopes = 1, plotflag = T
	)
{
#	This is a function which gives a plot of Higuchi's fractal
#	dimension algorithm for various
#	aggregation levels.  It calls the function Sfract2.S, which, in
#	turn, calls Cfract2.c.  The estimate of the slope, beta, gives
#	an estimate for the parameter, H, via :  H = 2 + beta
#	Variables:
#		data	--	The data being used.
#		nvar	-- 	The number of aggregation levels.
#		power1,power2-- Cut-offs for points to be used in estimating
#				beta.
#		slopes  --	Flag.  If 1, draw the estimated line whose 
#				slope is beta.
#		plotflag   --	Flag.  If T, plots the estimates vs. aggregation
#				level.
	cat("\n")	
	# n is the length of the data vector, absolute is the vector of length nvar+1
# which is the output of Sfract2.  It is the vector of curve lengths  at the 
# different levels of aggregation.
	n <- length(data)
	absolute <- Sfract2(data, nvar)	
	# The increment to be used in plotting and estimating the slope.  Depends on 
# the number of levels.
	increment <- (log10(n/2))/nvar	
	# lv and lm are going to be the y and x values respectively of the points in the
# plot.
	lv <- log10(absolute)
	lm <- rep(0, nvar)	
	# power1 and power2 determine whether the points will be included in the 
# estimate of beta.
	for(i in 1:nvar)
		lm[i] <- log10(floor(10^(i * increment)))	
	#Above line changed 2/28/95, to make the plotting consistent with calculations.
	lmadj <- lm[(lm >= power1) & (lm <= power2)]
	lvadj <- lv[(lm >= power1) & (lm <= power2)]
	lmadjc <- lm[(lm < power1) | (lm > power2)]
	lvadjc <- lv[(lm < power1) | (lm > power2)]	
	# The line is fitted to the appropriate points.
	fit <- lsfit(lmadj, lvadj)
	intercept <- fit$coef[1]
	beta <- fit$coef[2]	#Plot only if plotflag == T.  VT. 12/28/94.
	if(plotflag == T) {
		plot(c(0, lm), c(0, lv - lv[1]), xlab = "log10(m)", ylab = 
			"log10(Curve Length)", type = "n")
		points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.06)
		points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06)
		abline(0, -1.5, lty = 2)	
	#The least squares fitted line is drawn if slopes = 1.
		if(slopes == 1)
			lines(lmadj, intercept - lv[1] + (beta) * lmadj)
	}
#Output results to screen.
	cat(" beta = ", beta, "\n\n")
	cat("H = ", 2 + beta, "\n\n")	
	#returning a value for the function.VT. 12/28/94
	result <- beta + 2
	names(result) <- NULL
	return(result)
}


Sfract2 <- function(data, nvar)
{
	n <- length(data)
	if(!is.loaded(C.symbol("Cfract2")))
		dyn.load("Cfract2.o")
	.C("Cfractal",
		data = as.double(data),
		n = n,
		nvar = as.integer(nvar),
		v = double(nvar + 1))$v[2:(nvar + 1)]
}
