change1<-function(x1,a=1,b=5,zz=500,maxx=1860, pp=1, inter=500, boots=0,p=.95)
{
# This routine is for doing Intervention/Changepoint Analysis on unevenly sampled time series
# and it starts with a user-defined changepoint (zzz) and then calculates median of points that are within a fixed
# plus/minus distance (inter) of the changepoint (zzz). It then assesses the significance of the difference in the median by bootstrapping
# (boots) the data 1000 times (does this for each group independently) and assessing whether observed difference is outside of confidence interval (p=.95). 
# It will calculate change on absolute (pp=0) or percent basis (pp=1, difference as percent of median1)
#
# maxx = 1860. Was a cutoff I used to exclude data that post-Euroamerican settlement.
#
# STEP 1:
# Read in data and options
# txa = times for a series, txb=observations for a series
# Forces dataset to be a matrix and then takes reads in 
	dataa <- as.matrix(x1)
	xname1_name.cols(x1)
	txa_as.vector(dataa[,a])
	txb_as.vector(dataa[,b])
# STEP 2:  
# Counts the number of points +- zz and then sets up array in anticipation of bootstrap.
# inter1 is the lower (older if on AD scale) and inter1 is upper (younger if on AD scale date)
inter1_zz-inter
inter2_zz+inter
nrows_0
nrows_length(txa)
a1_0
a2_0
for (i in 1:nrows) {
   if(txa[i]<maxx){
        if(txa[i] <= zz & txa[i] >= inter1) {a1_a1+1}
        if(txa[i] <= inter2 & txa[i] >= zz) {a2_a2+1}
}
}

intera1_matrix(NA,nrow=a1,ncol=2)
intera2_matrix(NA,nrow=a2,ncol=2)
# STEP 3:
# Populate the intera1 and intera 2 matrices with original data values.  Column 1 is date, Column 2 is original values
	b1_0
	b2_0
for (i in 1:nrows) {
   if(txa[i]<maxx){
   if(txa[i] <= zz & txa[i] >= inter1) {
	b1_b1+1
	intera1[b1,1]_txa[i]
	intera1[b1,2]_txb[i]
                                       }
  if(txa[i] <= inter2 & txa[i] >= zz) {
	b2_b2+1
	intera2[b2,1]_txa[i]
	intera2[b2,2]_txb[i]
                                       }
					}
                    }

#STEP 4:
# Calculate the median difference between the two groups and save.  Then bootstrap values.  Note that there is a test below that checks if too
# many permutations are called for.
colss_boots+1
minboot_min(b1,b2)
minboot_minboot^minboot
if(colss > minboot) {colss_minboot+1}
outp_matrix(NA,1,colss)
time1_median(intera1[,1],na.rm=T)
time2_median(intera2[,1],na.rm=T)
median01_median(intera1[,2],na.rm=T)
median02_median(intera2[,2],na.rm=T)
tintera1_intera1[,2]
tintera2_intera2[,2]
for (j in 1:colss){
median1_median(tintera1,na.rm=T)
median2_median(tintera2,na.rm=T)
differ_median1-median2
if (pp == 1) {differ_(differ/median1)*100}
outp[1,j]_differ
tintera1_as.vector(sample(intera1[,2],replace=T))
tintera2_as.vector(sample(intera2[,2],replace=T))


}
# Collate/Collect Data and Save as Output Matrix (ciout)
#   ciout[,1]_Column of observations
#   ciout[,2]_change point
#   ciout[,3]_median time(txa) in group 1
#   ciout[,4]_median time(txa) in group 2
#   ciout[,5]_median of points in group 1
#   ciout[,6]_median of points in group 2
#   ciout[,7]_median1 - median2
#   ciout[,8] <- Lower CI for bootstrapped differences 
#   ciout[,9] <- Upper CI for bootstrapped differences
#   ciout[,10]_Number of points in group 1
#   ciout[,11]_Number of points in group 2

ciout <- matrix(NA, nrow = 1, ncol = 11)
   ciout[,1]_b   
   ciout[,2]_zz
   ciout[,3]_round(time1,0)
   ciout[,4]_round(time2,0)
   ciout[,5]_round(median01,3)
   ciout[,6]_round(median02,3)
   ciout[,7]_round(outp[1,1],4)
   ciout[,10]_a1
   ciout[,11]_a2
   ciout[,8] <- round(quantile(outp, probs=((1-p)/2), na.rm = T),4) 
   ciout[,9] <- round(quantile(outp, probs=(p+((1-p)/2)),na.rm=T),4)
ciout

# Close the Script
}