"""
Contains some functions for calculating cross-correlations, autocorrelations,
and Pearson correlations among rate units in a rate network
"""

using Statistics
import DSP.xcorr


"""
	xcorr_unbiased(x::Vector, y::Vector, maxlag::Int64)

Returns the cross-correlation of two signals normalized by the length of the
available convolution window.

Written by Rainer Engelken (re2365@columbia.edu).
"""
function xcorr_unbiased(x, y, maxlag)
	l = length(x) # we assume same length for x and y
	lags = -maxlag:maxlag
	scale = l .- abs.(lags)
	scale[scale .<= 0] .= 1 # avoid zero division
    xcorr(x, y)[l .+ lags] ./ scale
end



"""
	avg_autocorrs(A::Matrix{T}, maxlag::Int64, wait::Int64=0) where T <: AbstractFloat

Returns the average autocorrelation for the activity of each unit described in
the `N x NT` matrix `A`.

Note the autocorrelations are not mean-subtracted. The `maxlag` parameter
gives the maximum shift of the signal, and the `wait` parameter gives the number
of timesteps in the beginning to disregard.
"""
function avg_autocorrs(A::Matrix{T}, maxlag::Int64, wait::Int64=0) where T <: AbstractFloat

	(N, NT) = size(A)

	unit_autocorrs = mapslices(view(A, :, (wait + 1):NT), dims=2) do x
		return xcorr_unbiased(x, x, maxlag)
	end

	return vec(mean(unit_autocorrs, dims=1))
end
