# This routine is for doing Discrete Correlation Function on unevenly sampled time series
# and it is called from the rundcf function

# STEP 1:
# Read in data and options
# txa = times for a series, tyta=transformed values for a, txb=times for b series, tytb= transformed values for b.
# bin = minimum separation in time to calculate r for and 
# mlag equals max lag that statistic calculated for 
# nrun=1000 (default) is number of randomizations and p =0.95 for confidence interval 
# the total length of time series

dcf5 <-function(BP,txa,tyta,txb,tytb, bin,mlag,nrun,p,corr) {
options(object.size=1000000000)

# STEP 2:  
# Sets up arrays and dummy variables in anticipation of bootstrap.
		txa_as.matrix(txa)
		txb_as.matrix(txb)
		tya_as.matrix(tyta)
		tyb_as.matrix(tytb)

# Sets the bin size to median difference for the two series unless set by user
if (bin < 1)
	{
bin1_median(diff(txa[,1]))
bin2_median(diff(txb[,1]))
bin_ceiling(max(bin1,bin2))
	}
# Calculate number of discrete lag categories to calculate r for

nlags_1+(2*(ceiling(mlag/bin)))

# Adds the number of runs + 2 to provide for original run plus nrun bootstraps and then
# also dimension the rab matrix so can start with k=2 for first set of values. 
nruns_nrun+2

if (nruns > 1) {rab_matrix(NA,nrow=nlags,ncol=nruns)}

#Provide Labels for the midpoint of lag categories for rab matrix
label_0
c_1+((nlags-1)/2)

# First Loop here does label for lag = 0 out to max positive lag

for (i in 1:c) 	{
		rab[i,1]_label
	label_label+bin
					}
# Second Loop here does labels for lag starting with first to max neg lag
ca_c+1
label_0
for (i in ca:nlags) {
	label_label-bin
			rab[i,1]_label
	
						}

#Create the dab matrix, set max number of rows to number of pairs that have a lag within the set minlag and maxlag
#Boundaries specified/calculate above. Based on Years BP, so if proxy B lags A the event in B happens more recently which means A-B is positive.
#Dmax is the maximum distance and dmin is minimum distance that cross-correlation will be calculated for. 

pairsa_length(txa)
pairsb_length(txb)

dmax_max(rab[,1])+(bin/2)
dmin_dmax*-1
pair_0
	for(i in 1:pairsa)	{
	for(j in 1:pairsb)	{
		dabb_txa[i]-txb[j]
		if (dabb <= dmax & dabb > dmin) {
		pair=pair+1
	}}}
dab_ matrix(NA,nrow=pair,ncol=3)

# This starts main loop
# STEP 3a:
# Calculate all possible paired differences between two series where the difference in time is greater than min lag and greater than maxlag
# pair=0 is counter to increment which row is pairwise difference is being calculated and should have maxvalue=rows

for (k in 3:nruns) {
if (k < 4) {
pair=0
if(BP==1) {
	for(i in 1:pairsa)	{
	for(j in 1:pairsb)	{
		dabb_txa[i]-txb[j]
		if (dabb <= dmax & dabb > dmin) {
		pair=pair+1
		dab[pair,1]_dabb
		dab[pair,2]_tya[i]
		dab[pair,3]_tyb[j]
  											 }
							}
							}	
}

#This is identical to above loop but reverses subtraction to account for years AD vs BP
else if(BP==0) {
	for(i in 1:pairsa)	{
	for(j in 1:pairsb)	{
		dabb_txb[j]-txa[i]
		if (dabb <= dmax & dabb > dmin) {
		pair=pair+1
		dab[pair,1]_dabb
		dab[pair,2]_tya[i]
		dab[pair,3]_tyb[j]
  											 }
							}
							}	
}
}
# Step 3b:
# Now calculate r for the specified bin size and bins out to maxlag (positive and then negative).
correlate_matrix(NA,nrow=pair, ncol=2)
 for (i in 1:nlags) 	{
	cc_0
	winmin_rab[i,1]-(bin/2)
	winmax_rab[i,1]+(bin/2)
	for (j in 1:pair)			{
	dabb_dab[j,1]
	if (dabb > winmin & dabb <=winmax) 	{
		cc=cc+1
		correlate[cc,1]_dab[j,2]
		correlate[cc,2]_dab[j,3]
												}
								}
if (k==3){rab[i,2]_cc}
if(corr==0) {rab[i,k]_cor(correlate[1:cc,1],correlate[1:cc,2],na.method="omit")}
if(corr==1) {rab[i,k]_cor(rank(correlate[1:cc,1]),rank(correlate[1:cc,2]),na.method="omit")}

						}
#STEP 4:
if (nrun> 1) {
	dab[,2]_as.matrix(sample(as.matrix(tya),size=length(dab[,2]),replace=T))
	dab[,3]_as.matrix(sample(as.matrix(tyb),size=length(dab[,3]),replace=T))
				}
# Close Loop for Bootstrap (k)
	}

# Step 5:
# Assess the statistical significance of the observed bootstrap values and then prepare a matrix with observed r values
# and upper and lower confidence interval.
ciout <- matrix(NA, nrow = nrow(rab), ncol = 5)
   ciout[,1]_rab[,1]
   ciout[,2]_rab[,3]
   ciout[,5]_rab[,2]
nrws=nrow(rab)
for(i in 1:nrws) {
#	Note that the -1 tells program to skip first column but do quantile for all other columns
    	ciout[i,3] <- quantile(rab[i,-(1:2)], probs=((1-p)/2), na.rm = T) 
       ciout[i,4] <- quantile(rab[i,-(1:2)], probs=(p+((1-p)/2)),na.rm=T)
				
}
#print(memory.size())
rab_0
correlate_0
dabb_0
tya_0
tyb_0
#print(memory.size())	
res <- list(bin=ciout[,1],rbin = ciout[,2],lower=ciout[,3], upper=ciout[,4],numbin=ciout[,5])
#Closes the Script
}