pro varadi_csa,filex,filey,rmax,M,K,filters,seed
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; NAME:
;	varadi
; CALLING SEQUENCE:
;	varadi_csa,filex,filey,rmax,M,K,filters,seed
; PURPOSE:
;	Compute the RLSSA from Varadi et al, ApJ, 1999, 526, 1052
; INPUTS:
;	filex		filename for first series
;	filey		filename for second series
;	rmax		number of largest eigenvalues for reconstruction
;	M		number of random lags
;	K		maximum lag
;	filters		reconstructed filters of size fltarr(K)
;	seed		dummy argument
; OPTIONAL KEYWORDS:
; 	none
; OUTPUTS:
;	none
; OPTIONAL OUTPUT:
;	none
; EXAMPLE:
;	thomson,file
; LIMITATIONS:
;	None yet
; COMMONS:
;	None
; PROCEDURES USED:
;       fft, readfits
; MODIFICATION HISTORY:
;	Beta 1: WRITTEN, Thierry Appourchaux, November 24, 2000
;------------------------------------------------------------------------------
;
; Read files
	print,'r_max (10 is a good number)=',rmax
	print,'M (number of different lags, 250 is OK)',M
	print,'K (maximum lag, 25000 is OK)',K
	a_x=readfits(filex,/silent)
	m_2=moment(a_x,/double)
	a_x=(a_x-m_2(0))/sqrt(m_2(1))
	a_y=readfits(filey,/silent)
	m_2=moment(a_y,/double)
	a_y=(a_y-m_2(0))/sqrt(m_2(1))
	!p.multi=[0,1,2]
;
; Check if x and y are the same, if yes flag=1
;
	diff=total(a_x-a_y)
	flag=0
	if (diff eq 0.) then flag=1
	print,diff,flag

;	a=randomn(seed,131072)*1.d0
;
;
	N_a=N_elements(a_x)
;
	k_lag=K
;
; Here we assume that the time series is sampled at 60 s
;
	xfreq=findgen(N_a)*1.e06/60./N_a
	a_rev_x=reverse(a_x)
	a_rev_y=reverse(a_y)
;
;	Pad the time series with zero
;
	a_new_x=dblarr(2*N_a)
	a_new_rev_x=dblarr(2*N_a)
	a_new_x(0:N_a-1)=a_x
	a_new_rev_x(0:N_a-1)=a_rev_x
;
	a_new_y=dblarr(2*N_a)
	a_new_rev_y=dblarr(2*N_a)
	a_new_y(0:N_a-1)=a_y
	a_new_rev_y(0:N_a-1)=a_rev_y
;
;	Then compute the correlation using fft (3.3)
;	Please note the use of the different normalization
;	It was calibrated using the real formula (3.1)
;

	q=2.d0*N_a*fft(fft(a_new_x,-1,/double)*fft(a_new_rev_y,-1,/double),1,/double)
	q_x=2.d0*N_a*fft(fft(a_new_x,-1,/double)*fft(a_new_rev_x,-1,/double),1,/double)
	q_y=2.d0*N_a*fft(fft(a_new_y,-1,/double)*fft(a_new_rev_y,-1,/double),1,/double)
;
; Shift to put the max at 0
;
	q=shift(double(q),-(N_a-1))
	q_x=shift(double(q_x),-(N_a-1))
	q_y=shift(double(q_y),-(N_a-1))
	qt=dblarr(2L*k_lag-1)
	qt_x=dblarr(2L*k_lag-1)
	qt_y=dblarr(2L*k_lag-1)
	qt(0:k_lag-1)=q(0:k_lag-1)
	qt(k_lag:2L*k_lag-2)=q(2L*N_a-1-(k_lag-2):2L*N_a-1)
;
	qt_x(0:k_lag-1)=q_x(0:k_lag-1)
	qt_x(k_lag:2L*k_lag-2)=q_x(2L*N_a-1-(k_lag-2):2L*N_a-1)
;
	qt_y(0:k_lag-1)=q_y(0:k_lag-1)
	qt_y(k_lag:2L*k_lag-2)=q_y(2L*N_a-1-(k_lag-2):2L*N_a-1)
;
; Calibration
;
;
;	q_clas=fltarr(k_lag)
;	for k=0,k_lag-1 do begin
;		q_clas(k)=total(a(0:N_a-1-k)*a(k:N_a-1))
;	endfor
;
; Equation (3.5) corrected...
;
	norm=N_a-abs(findgen(2L*k_lag-1)-(k_lag-1))
	qt=qt/double(norm)
	qt_x=qt_x/double(norm)
	qt_y=qt_y/double(norm)
;
; normalize, equation (58)
;
	if (flag eq 1) then begin
		c=qt/qt(0)
		c_x=qt_x/qt_x(0)
		c_y=qt_y/qt_y(0)
	endif else begin
		c=qt
		c_x=qt_x
		c_y=qt_y
	endelse

;
; Set randomly the lag and makes sure that they are all different
;
	C_diff=dblarr(M,M)
	C_diff_x=dblarr(M,M)
	C_diff_y=dblarr(M,M)
	S_x=floor(k_lag*randomu(seed,M))+1
	S_y=floor(k_lag*randomu(seed,M))+1
	if (flag eq 1) then S_y=S_x
	for i=0,M-1 do begin
		for j=i+1,M-1 do begin
			;print,i,j
			C_diff_x(i,j)=S_x(i)-S_x(j)
			C_diff_y(i,j)=S_y(i)-S_y(j)
			repeat begin
				nodiff=1
				if ((C_diff_y(i,j) eq 0.) or (C_diff_x(i,j) eq 0.)) then begin
					print,C_diff_x(i,j),C_diff_y(i,j),i,j,S_x(i),S_x(j),S_y(i),S_y(j)
					S_x(j)=floor(k_lag*randomu(seed))+1
					S_y(j)=floor(k_lag*randomu(seed))+1
					if (flag eq 1) then S_y(j)=S_x(j)
					C_diff_x(i,j)=S_x(i)-S_x(j)
					C_diff_y(i,j)=S_y(i)-S_y(j)
					print,C_diff_x(i,j),C_diff_y(i,j),C_diff(i,j),i,j
					nodiff=0
				endif
			endrep until nodiff
		endfor
	endfor
	for i=0,M-1 do begin
		for j=0,M-1 do begin
			C_diff(i,j)=S_y(i)-S_x(j)
		endfor
	endfor

;
; Set the random matrix
;
	C_rand=dblarr(M,M)
	C_rand_x=dblarr(M,M)
	C_rand_y=dblarr(M,M)
	for i=0,M-1 do begin
		for j=i,M-1 do begin
			C_rand_x(i,j)=c_x((2L*k_lag-1)*(1-sign(C_diff_x(i,j)))*0.5+C_diff_x(i,j))
			C_rand_y(i,j)=c_y((2L*k_lag-1)*(1-sign(C_diff_y(i,j)))*0.5+C_diff_y(i,j))
			C_rand_x(j,i)=C_rand_x(i,j)
			C_rand_y(j,i)=C_rand_y(i,j)
		endfor
	endfor
	for i=0,M-1 do begin
		for j=0,M-1 do begin
			C_rand(i,j)=c((2L*k_lag-1)*(1-sign(C_diff(i,j)))*0.5+C_diff(i,j))
		endfor
	endfor
	Print,'Eigenvalues...'
;
; Compute SVD according to Eq (2.3)
;
	svdc,C_rand,lambda,evecs_x,evecs_y,/double
;
; Compute sigma_x according to (2.5)
;
	pseudo_x=transpose(evecs_x) ## C_rand_x ## evecs_x
	lambda_x=dblarr(M)
	for i=0,M-1 do begin
		lambda_x(i)=pseudo_x(i,i)
	endfor
;
; Compute sigma_x according to (2.5) modified for y
;
	lambda_y=dblarr(M)
	pseudo_y=transpose(evecs_y) ## C_rand_y ## evecs_y
	for i=0,M-1 do begin
		lambda_y(i)=pseudo_y(i,i)
	endfor
;
; Set the number of eigenvalues (or highest correlation coefficient)
;
	r=rmax
;
; Compute B according to (3.9)
;
	ratio=lambda/lambda_x
;
; Compute correlation coefficient according to (3.10)
;
	r_corr=lambda/sqrt(lambda_x)/sqrt(lambda_y)
;
; Sort data
;
	if (flag eq 1) then begin
		index_c=reverse(sort(lambda))
	endif else begin
		index_c=reverse(sort(r_corr))
	endelse
	evecs=1.d0
	eigenvalues=eigenql(C_rand_x,eigenvectors=evecs,/double)
;
; Set filter vector to 0
;
	f1=dblarr(k_lag)
;
; Compute filter
;
	evecs_x=transpose(evecs_x)
	evecs_y=transpose(evecs_y)
	Print,'Filter...'
	for mm=0,r-1 do begin
		for i=0,M-1 do begin
			nn=S_x(i)-S_y
			index=where(nn ge 0.)
			if (index(0) ne -1) then begin
			f1(nn(index))=f1(nn(index))+evecs_x(i,index_c(mm))*evecs_y(index,index_c(mm))*ratio(index_c(mm))/M
			endif
		endfor
	endfor
	f1(0)=0.
	filters=f1
end