bam.func <- function(sitesdat, radius, cliquesdat)
{
	# This is the main function for the BAM algorithm.
	#
	# DEMONSTRATION EXAMPLE: 
	#
	# 1. Create a pair of maps to be matched, each on a 40 x 40 grid of 
	# squares, with each square of the grid having presence with 
	# probability 0.2 (and absence with probability 0.8):
	# 	mymaps <- makemap.func(40, 40, 0.2)
	#
	# 2. Find the best attainable match of these maps with cliques as the
	# eight nearest neighbours:
	# (note: this means "radius" should be sqrt(2), where a small number is
	# added to avoid rounding errors):
	# 	bam.func(mymaps, radius=sqrt(2)+0.01)
	#
	# ALTERNATIVE: (for the same results: following from step 1 above):
	# 
	# 2. Create the cliques data manually for the eight nearest neighbours:
	# 	mycliquesdat <- makeclique.func(mymaps, radius=sqrt(2)+0.01)
	#
	# 3. Find the best attainable match:
	# 	bam.func(mymaps, cliquesdat=mycliquesdat)
	#
	# This alternative approach allows you to store and view the cliques
	# data, and it will make the program run a little faster if you need
	# to match several different maps on the same grid.  (The makeclique.func
	# routine only needs to be run once for a given grid set-up, even if you
	# generate completely new maps: 
	# for example, the following lines only require the makeclique line once:
	# 	mymaps1 <- makemap.func(40, 40, 0.5)
	#	mymaps2 <- makemap.func(40, 40, 0.1) # (must be a 40x40 grid again)
	# 	mycliquesdat <- makeclique.func(mymaps1, radius=1.5)
	# 	bam.func(mymaps1, cliquesdat=mycliquesdat)
	# 	bam.func(mymaps2, cliquesdat=mycliquesdat)
	# 
	# The results should be the same as
	# 	bam.func(mymaps1, radius=1.5)
	# 	bam.func(mymaps2, radius=1.5)
	# 
	# The makeclique option is also helpful to see what is required if you need 
	# to enter cliques manually (rather than using a radius idea).  To view the 
	# cliques data, just type  
	# 	mycliquesdat
	# See also the General Information below for more details.
	#
	# GENERAL INFORMATION:
	# 
	# sitesdat is a data frame with the following headers:
	# site, map1, map2 (the presence / absence in the first and second surveys), 
	# x, y (site coordinates)
        # There is a choice of entering clique data manually, or of using a radius
        # option so that each site's clique contains all sites within the given 
	# radius. Only one of the two arguments "radius" and "cliquesdat" should
	# be provided, otherwise the program will stop. 
	# Using radius = sqrt(2) (plus a small number to avoid rounding errors)
	# gives an eight-nearest-neighbour result for sites on a regular grid with 
	# sites placed exactly one unit apart.
        #
        # If the cliques are to be read in from dataframe cliquesdat, the clique 
	# data should be arranged as follows.  Suppose that
        # site 1 has clique (2, 3, 4);
        # site 2 has clique (1, 10);
        # site 3 has clique (1, 4, 21); and so on.
        # Then cliquesdat should be a data frame as so:
        #
        # site clique
        # 1    2
        # 1    3
        # 1    4
        # 2    1
        # 2    10
        # 3    1
        # 3    4
        # 3    21
        #
        # etc.
	# 
	# First some housekeeping: only one of cliquesdat and radius should  
        # be provided:
        if((!missing(cliquesdat))&(!missing(radius)))
                stop("Only one of radius and cliquesdat should be provided.")
        # 
	map1 <- sitesdat$map1
	map2 <- sitesdat$map2
	xvec <- sitesdat$x
	yvec <- sitesdat$y
	Nsites <- length(map1)
	mismatch0 <- (1:Nsites)[(map1 == 0) & (map2 == 1)]
	mismatch1 <- (1:Nsites)[(map1 == 1) & (map2 == 0)]
	mismatchall <- sort(c(mismatch0, mismatch1))
	Nmism0 <- length(mismatch0)
	Nmism1 <- length(mismatch1)
	Nmism <- Nmism0 + Nmism1
	if(!missing(cliquesdat))	
		swapmat <- clique.func(mismatchall, mismatch0, mismatch1, sitesdat,
						cliquesdat=cliquesdat)
	else swapmat <- clique.func(mismatchall, mismatch0, mismatch1, sitesdat,
		radius = radius)
	swapsing <- 0
	n2blocks <- 0
	n3blocks <- 0
	n4blocks <- 0
	n5blocks <- 0
	nblocks <- 0
	swapperf <- 0
	sing1 <- singletons.func(swapmat)
	swap <- sing1$swap
	swapmat <- data.frame(sing1$swapmat)
	swapsing <- sing1$swap
	catline("Singletons sweep 1 completed.  Swap =", swapsing)
	# If swapmat has 0 rows, the BAM procedure is finished.
	# If not, go into the connected clique, perfect cluster, and
	# block-swapping routines.
	# Note that "perfect clusters" and "exact connected components"
	# are the same thing (just different notational eras!)
	while(nrow(swapmat) != 0) {
		perfect <- T
		while((nrow(swapmat) != 0) & perfect) {
			perfect <- F
			selsite <- min(swapmat$site)
			adjmat <- connclique.func(swapmat, selsite)
			# catline("Selected site: ", selsite)
			# print(adjmat)
			nrow.adjmat <- nrow(adjmat)
			if(nrow.adjmat == ncol(adjmat)) {				
				perfect <- checkperfect.func(adjmat)
				# if adjmat was perfect it can be swapped
				# and all sites in the connected clique must
				# be deleted.
				# full swap is achieved for a perfect cluster.
				# full swap = nrow.adjmat = ncol.adjmat
				if(perfect) {
					swap <- swap + nrow.adjmat
					swapperf <- swapperf + nrow.adjmat
					swapsites <- as.numeric(unlist(dimnames(
						adjmat)))
					swapmat <- swapmat[!is.element(swapmat$
						site, swapsites),  ]
					catline("Exact connected component: swap =", 
						nrow.adjmat)
				}
			}
			# catline("Searching for perfect clusters: ", perfect)
		}
		if(nrow(swapmat) != 0) {
			# begin searching for blocks and performing singletons 
			# search for two-blocks in adjmat
			twob <- twoblock.func(adjmat)
			if(twob$status == "none")
				twob <- twoblock.func(t(adjmat))
			# catline("Searching for 2-blocks: ")
			# print(twob)
			if(twob$status == "found") {
				# swap the block: swap=2
				swap <- swap + 2
				# delete the sites in the block from swapmat, 
				# and also remove them from the "swap" set of
				# any other sites in swapmat:
				swapsites <- as.numeric(unlist(twob$block))
				swapmat <- swapmat[!is.element(swapmat$site,
					swapsites),  ]
				swapmat <- swapmat[!is.element(swapmat$swap,
					swapsites),  ]
				n2blocks <- n2blocks + 1
				nblocks <- nblocks + 1
				catline("Swapped 2-block")
			}
			else {
				r <- 2
				foundblock <- F
				while(!foundblock) {
					r <- r + 1
					# search for (r, r)-blocks in adjmat
					rsb <- rsblock.wrap(adjmat, r)
					if(rsb$status == "none")
						rsb <- rsblock.wrap(t(adjmat),
							r)
					# catline("Searching for ", r, "-blocks: ")
					# print(rsb)
					if(rsb$status == "found")
						foundblock <- T
				}
				# the above "while" loop exits when a block is found:
				# the block is stored in rsb$block
				# swap the block: swap=r
				catline("Swapped ", r, "-block")
				swap <- swap + r
				# delete the sites in the block from swapmat, 
				# and also remove them from the "swap" set of
				# any other sites in swapmat:
				swapsites <- as.numeric(unlist(rsb$block))
				swapmat <- swapmat[!is.element(swapmat$site,
					swapsites),  ]
				swapmat <- swapmat[!is.element(swapmat$swap,
					swapsites),  ]
				nblocks <- nblocks + 1
				if(r == 3)
					n3blocks <- n3blocks + 1
				if(r == 4)
					n4blocks <- n4blocks + 1
				if(r == 5)
					n5blocks <- n5blocks + 1
			}
		}
		if(nrow(swapmat) != 0) {
			# do singletons again
			singnew <- singletons.func(swapmat)
			swap <- swap + singnew$swap
			swapmat <- data.frame(singnew$swapmat)
			if(singnew$swap>0) 
				catline("Singletons new sweep completed.  Swap =", singnew$swap)
			swapsing <- swapsing + singnew$swap
		}
	}
	BAM <- Nsites - Nmism + 2 * swap
	catline()
	catline()
	catline("SUMMARY")
	catline("Swap attributable to singletons: ", swapsing)
	catline("Swap attributable to exact connected components: ", swapperf)
	catline("Number of 2-blocks: ", n2blocks)
	catline("Number of 3-blocks: ", n3blocks)
	catline("Number of 4-blocks: ", n4blocks)
	catline("Number of 5-blocks: ", n5blocks)
	catline("Total number of blocks: ", nblocks)
	catline()
	catline("Number of sites: ", Nsites)
	catline("Original match: ", Nsites - Nmism)
	catline("Best possible swap: ", swap)
	catline("Best attainable match: ", BAM)
	catline()
	catline()
	BAM
}
checkperfect.func <- function(adjmat)
{
	# checkperfect.func
	# checks whether a connected clique is a perfect cluster.
	valency <- unique(c(apply(adjmat, 1, sum), apply(adjmat, 2, sum)))
	perfect <- length(valency) == 1
	# returns logical flag: perfect or not
	perfect
}
clique.func <- function(mismsites, mism0sites, mism1sites, sitesdat, radius, 
cliquesdat)
{
	# clique.func
	# Takes the list of mismatched site numbers, mismsites, and constructs a
	# data frame containing all mismatches and swappable sites (mismatched sites
	# with the opposite status in the same clique).
	# 
	# There is a choice of entering clique data manually, or of using a radius
	# option so that each site's clique contains all sites within the given radius.
	# The default is to use radius = sqrt(2) (plus a small number to avoid rounding 
	# errors).  This gives an eight-nearest-neighbour result for sites on a regular
	# grid with sites placed exactly one unit apart.
	# 
	# If the cliques are to be read in from dataframe cliquesdat, the clique data 
	# should be arranged as follows.  Suppose that 
	# site 1 has clique (2, 3, 4);  
	# site 2 has clique (1, 10);
	# site 3 has clique (1, 4, 21); and so on.  
	# Then cliquesdat should be a data frame as so: 
	#
	# site clique
	# 1    2
	# 1    3
	# 1    4
	# 2    1
	# 2    10
	# 3    1
	# 3    4
	# 3    21
	#
	# etc.  
	# The output from this function (clique.func) is a new data frame with only the
	# mismatched sites listed in the first column, and only the swappable 
	# clique-sites listed in the second.  For example (using the cliques above): if 
	# site 1 is a mismatched 0; 
	# site 2 is a match; 
	# site 3 is a mismatched 0; 
	# site 4 is a mismatched 1; 
	# site 21 is a mismatched 1;
	# then the output would be:
	#
	# site swap
	# 1    4
	# 3    4
	# 3    21
	# etc.
	#
	#
	# If there are no swappable sites, the mismatched site does not appear in
	# the output.  If the site is a match, it does not appear in the output either.
	#
	# first some housekeeping: only one of cliquesdat and radius should
	# be provided:
	if((!missing(cliquesdat))&(!missing(radius)))
		stop("Only one of radius and cliquesdat should be provided.")
	#	
	# now introduce the function for finding clique sites for site i
	# using the radius option
	#
	nsites <- length(sitesdat$map1)
	mismdat.sites <- (1:nsites)[is.element((1:nsites), mismsites)]
	mismdat <- sitesdat[mismdat.sites,]
	mismdat$site <- mismdat.sites
	cliqueset.func <- function(i)
	{
		# print(i)
		i.x <- mismdat$x[mismdat$site == i]
		i.y <- mismdat$y[mismdat$site == i]
		dist <- sqrt((mismdat$x - i.x)^2 + (mismdat$y - i.y)^2)
		clique <- mismdat$site[dist <= radius]
		n.clique <- length(clique)
		as.vector(rbind(rep(i, n.clique), clique))
	}
	if(!missing(cliquesdat)) {
		# Transform the cliquesdat data frame so that sites are 
		# numbered 1:nsites, even if the user's numbering was
		# not originally 1:nsites:
		# the result orders site 1 to site "nsites" in the order
		# that they appear in sitesdat$site.
		cliquesdat$site <- match(cliquesdat$site, sitesdat$site)	
		cliquesdat$clique <- match(cliquesdat$clique, sitesdat$site)	
		# now continue with transformed cliquesdat:
		len <- length(cliquesdat$site)
		# extract the "site" rows involving mismatched sites
		output <- cliquesdat[(1:len)[is.element(cliquesdat$site, 
			mismsites)],  ]
		# extract the "clique" rows involving mismatched sites
		len2 <- length(output$site)
		output <- output[(1:len2)[!is.na(match(output$clique, mismsites
			))],  ]
		# remove the rows with a mism0 site and a mism0 clique
		len3 <- length(output$site)
		output <- output[(1:len3)[!(is.element(output$site, mism0sites) &
			is.element(output$clique, mism0sites))],  ]
		# remove the rows with a mism1 site and a mism1 clique
		len4 <- length(output$site)
		output <- output[(1:len4)[!(is.element(output$site, mism1sites) &
			is.element(output$clique, mism1sites))],  ]
		names(output) <- c("site", "swap")
		# add a column to say whether the site is a mism0 or mism1
		len5 <- length(output$site)
		output$type <- rep(0, len5)
		output$type[is.element(output$site, mism1sites)] <- 1
	}
	else {
		# this is the case where cliquesdat is missing so we use the 
		# radius option
		output <- matrix(unlist(sapply(mismsites, cliqueset.func)),
			ncol = 2, byrow = T)
		output <- data.frame(site = output[, 1], clique = output[,
			2])
		# the next lines are identical to the cliquesdat case.
		# remove the rows with a mism0 site and a mism0 clique
		len3 <- length(output$site)
		output <- output[(1:len3)[!(is.element(output$site, mism0sites) &
			is.element(output$clique, mism0sites))],  ]
		# remove the rows with a mism1 site and a mism1 clique
		len4 <- length(output$site)
		output <- output[(1:len4)[!(is.element(output$site, mism1sites) &
			is.element(output$clique, mism1sites))],  ]
		names(output) <- c("site", "swap")
		# add a column to say whether the site is a mism0 or mism1       
		len5 <- length(output$site)
		output$type <- rep(0, len5)
		output$type[is.element(output$site, mism1sites)] <- 1
	}
	output
}
connclique.func <- function(swapmat, selsite)
{
	# connclique.func
	# Finds the connected clique of the selected site, selsite.
	# Returns the adjacency matrix of the connected clique, with 
	# site numbers as the row and column names, and 0's in the rows, 
	# 1's in the columns. 
	site <- swapmat$site
	swap <- swapmat$swap
	type <- swapmat$type
	clique <- selsite
	clique <- c(selsite, swap[site == selsite])
	addsites <- unique(swap[is.element(site, clique)])
	while(any(!is.element(addsites, clique))) {
		clique <- unique(c(clique, addsites))
		addsites <- unique(swap[is.element(site, clique)])
	}
	clique <- sort(clique)
	# now get into adjacency matrix format
	clique.type <- type[match(clique, site)]
	# clique0 and clique1 are the site numbers of the 0 and 1 sites
	# in the connected clique.
	clique0 <- clique[clique.type == 0]
	clique1 <- clique[clique.type == 1]
	nrows <- length(clique0)
	ncols <- length(clique1)
	adjmat <- matrix(0, nrow = nrows, ncol = ncols)
	dimnames(adjmat) <- list(clique0, clique1)
	for(i in 1:length(clique1)) {
		s <- clique1[i]
		adjmat[, i] <- as.numeric(is.element(clique0, swap[site == s]))
	}
	adjmat

}
makeclique.func <- function(mapsdat, radius)
{
	# makeclique.func 
	# forms a "cliquesdat" data file based on radius, suitable for reading 
	# into the bam.func program.  
	# EXAMPLE: 
	# 	mymaps <- makemap.func(20, 20)
	#	mycliques <- makeclique.func(mymaps, radius=2)
	# 	bam.func(mymaps, cliquesdat=mycliques)
	# This has the same effect as:
	# 	bam.func(mymaps, radius=2)
	# but with the additional effect that the clique structure is then
	# stored in object mycliques. 
	#
	inner.func <- function(s){
		rbind(sites[s], mapsdat$site[(mapsdat$x-mapsdat$x[s])^2 + 
		(mapsdat$y-mapsdat$y[s])^2 <= radius^2])
	}
	sites <- mapsdat$site
	nsites <- length(sites)
	temp <- sapply(1:nsites, inner.func)
	cliquesdat <- matrix(unlist(temp), ncol=2, byrow=T)
	cliquesdat <- cliquesdat[cliquesdat[,1]!=cliquesdat[,2], ]
	cliquesdat <- as.data.frame(cliquesdat)
	names(cliquesdat) <- c("site", "clique")
	cliquesdat
}
makemap.func <- function(nx, ny, p.pres = 0.5)
{
	# makemap.func
	# constructs a pair of presence / absence distributions on a grid of 
	# dimensions nx by ny, with P(presence) = p.pres for every site.
	grid <- expand.grid(x = 1:nx, y = 1:ny)
	nsites <- nx * ny
	grid$site <- 1:nsites
	grid$map1 <- rbinom(nsites, 1, p.pres)
	grid$map2 <- rbinom(nsites, 1, p.pres)
	grid
}
rsblock.func <- function(adjmat, r, s, rsites, ssites, adjcolsites, adjrowsites, 
subb)
{
	# rsblock.func
	# searches in adjmat (the adjacency matrix) for a block of r columns which 
	# account for the total valency of *at least* s rows.
	# (Note that the notation is slightly different from the Biometrics paper:
	# there, an (r, s)-block had r columns and exactly s rows; here, it has
	# r columns and at least s rows.)
	# rsites and ssites are the SITE NUMBERS of the sites currently in the
	# block, at the time of the recursive call.
	# adjrowsites and adjcolsites are numeric vectors of the SITE NUMBERS of 
	# sites in the matrix adjmat: so adjrowsites[i] is the site number of the 
	# site corresponding to row i in adjmat.
	# subb stands for "subblocks": subb=F if it is known that there are
	# no subblocks in adjmat; otherwise subb=T.
	#
	# Note: the adjacency matrix does not operate with row and column
	# names (unlike the rsblock.verbose function) for efficiency.
	# The site numbers attached to the rows and columns are kept in 
	# adjrowsites and adjcolsites instead.
	#
	status <- "cont"
	# ("cont" for continue)
	if(length(adjrowsites) < s) status <- "none"
	if(length(adjcolsites) < r)
		status <- "none"
	if(status == "cont") {
		# find valencies of rows:
		if(length(adjcolsites) == 0) {
			rowval <- rep(0, length(adjrowsites))
		}
		else {
			adjmat <- matrix(adjmat, ncol = length(adjcolsites))
			rowval <- apply(adjmat, 1, sum)
		}
		# delete rows with valency > r : leaves nrows rows and ncols columns
		adjrowsites <- adjrowsites[rowval <= r]
		nrows <- length(adjrowsites)
		ncols <- length(adjcolsites)
		if(length(adjmat) > 0) {
			adjmat <- matrix(adjmat[rowval <= r,  ], nrow = nrows,
				ncol = ncols)
		}
		if(length(adjrowsites) < s)
			status <- "none"
	}
	if((status == "cont") & (r == 0)) {
		# base case where r=0:
		# if nrows<s then there is no block, 
		# otherwise a block is found and is now constructed
		if(nrows < s) status <- "none" else {
			status <- "found"
			ssites <- c(ssites, adjrowsites)
			block <- list(r = rsites, s = ssites)
		}
	}
	if((status == "cont") & (r > 0)) {
		# case r > 0 
		if(length(adjrowsites) == 0) colval <- rep(0, length(
				adjcolsites)) else colval <- apply(adjmat,
				2, sum)
		if(subb == F) {
			# delete columns with 1 entry, as they can't be part
			# of a block with no subblocks.
			# Also delete the associated row.
			ncols <- length(adjcolsites)
			nrows <- length(adjrowsites)
			deletecols <- (1:ncols)[colval == 1]
			# quick way of finding rows to delete:
			if(length(deletecols) > 0) {
				deleterows <- matrix(1:nrows, nrow = 1) %*%
					adjmat[, deletecols]
				deleterows <- unique(as.vector(deleterows))
				ncols <- ncols - length(deletecols)
				adjmat <- matrix(adjmat[,  - deletecols], ncol
					 = ncols, nrow = nrows)
				adjcolsites <- adjcolsites[ - deletecols]
				colval <- colval[ - deletecols]
				nrows <- nrows - length(deleterows)
				adjmat <- matrix(adjmat[ - deleterows,  ],
					ncol = ncols, nrow = nrows)
				adjrowsites <- adjrowsites[ - deleterows]
			}
		}
		# now delete columns with no entries (for subb=T or F)
		adjmat <- as.matrix(adjmat[, colval > 0])
		adjcolsites <- adjcolsites[colval > 0]
		if(length(adjcolsites) < r)
			status <- "none"
		if(length(adjrowsites) < s)
			status <- "none"
	}
	if(status == "cont") {
		# status="cont" at this stage implies that r>0.
		# Find selected row with maximum valency to add to block:
		# selrowpos is the row position in adjmat;
		# selrowsite is the site number of the selected row.
		rowval <- apply(adjmat, 1, sum)
		nrows <- length(adjrowsites)
		ncols <- length(adjcolsites)
		selrowpos <- min((1:nrows)[rowval == max(rowval)])
		selrowsite <- adjrowsites[selrowpos]
		selrowval <- rowval[selrowpos]
		# selrowval is the valency of the selected site.
		# Next add selrowsite and its intersecting columns to the block:
		newssites <- c(ssites, selrowsite)
		selcolpos <- (1:ncols)[adjmat[selrowpos,  ] == 1]
		selcolsites <- adjcolsites[selcolpos]
		newrsites <- c(rsites, selcolsites)
		# set up new matrix with the selected row and columns deleted:
		newrowsites <- adjrowsites[ - selrowpos]
		newcolsites <- adjcolsites[ - selcolpos]
		newmat <- matrix(adjmat[ - selrowpos,  ], ncol = ncols)
		newmat <- matrix(newmat[,  - selcolpos])
		blockcall <- rsblock.func(newmat, r - selrowval, s - 1, 
			newrsites, newssites, newcolsites, newrowsites, subb = 
			T)
		# if there is no block containing selrow then delete selrow
		# from adjmat
		if(blockcall$status == "none") {
			adjmat <- matrix(adjmat[ - selrowpos,  ])
			adjrowsites <- adjrowsites[ - selrowpos]
			adjmat <- matrix(adjmat, ncol = length(adjcolsites))
			blockretry <- rsblock.func(adjmat, r, s, rsites, ssites,
				adjcolsites, adjrowsites, subb = F)
			block <- blockretry$block
			status <- blockretry$status
		}
		else {
			block <- blockcall$block
			status <- blockcall$status
		}
	}
	if(status == "none")
		block <- NULL
	list(block = block, status = status)
}
singletons.func <- function(swapmat)
{
	# singletons.func
	#
	# There are two types of site numbering: 
	# firstly the site number, which is the absolute site number 
	# decided before the maps are observed;
	# secondly, the list number, which is the position of the site
	# in the matrix swapmat. 
	# For example, if swapmat is as follows:
	#  site swap type
	#     2    6    0
	#     2    8    0
	#     3    8    0
	#     5    2    1
	#     5   12    1
	#     8    2    1
	#     8    4    1
	# then the list numbers would be
	#   list number    site number 
	#          1		  2
	# 	   2 		  3
	#          3              5
	# 	   4              8
	# etc.
	# (The list numbers are always 1, 2, 3, 4, ..., n.sites.  
	# To transfer from the list number to site number and back:
	# siteno = site.number[listno]
	# listno = (1:n.sites)[site.number==siteno]
	# or (for a vector of siteno's):
	# listno = (1:n.sites)[match(sitenos, site.number)] 
	n.entries <- length(swapmat$site)
	site.index <- (1:n.entries)[!duplicated(swapmat$site)]
	site.number <- swapmat$site[site.index]
	# (same as unique(swapmat$site))
	site.type <- swapmat$type[site.index]
	site.val <- tapply(swapmat$swap, swapmat$site, length)
	n.sites <- length(site.number)
	# site.index[i] is the row number in matrix swapmat for the first 
	# appearence of list site i.
	# site.number[i] is the site number of list site i;
	# site.type[i] is the type (mismatched 0 or 1) of list site i;
	# site.val[i] is the valency of list site i.
	swap.listno <- NA
	if(any(site.val==1)) swap.listno <- min((1:n.sites)[site.val == 1])
	# swap.listno is the list number of the site to be swapped;
	# swap.siteno is the site number of the site to be swapped.
	n.swaps <- 0
	while(!is.na(swap.listno)) {
		swap.siteno <- site.number[swap.listno]
		swapwith.siteno <- swapmat$swap[swapmat$site == swap.siteno]
		swapwith.listno <- (1:n.sites)[site.number == swapwith.siteno]
		# swap.siteno is the site number of the singleton;
		# swapwith.siteno is the site number of the single site 	
		# it can be swapped with.
		# swapwith.listno is the list number of the single site it can
		# be swapped with.
		# 
		# set the valency of the swapwith site to zero, then reduce by 
		# 1 the valency of all sites valent with swapwith (this 
		# includes the swap site);
		# the sites valent with swapwith are called knockon.
		site.val[swapwith.listno] <- 0
		knockon.siteno <- swapmat$swap[swapmat$site == swapwith.siteno]
		knockon.listno <- (1:n.sites)[match(knockon.siteno, site.number
			)]
		site.val[knockon.listno] <- site.val[knockon.listno] - 1
		
		# remove all records of sites valent with swapwith.siteno
		# (this includes removing the original site, swap.siteno):
		swapmat <- swapmat[swapmat$swap != swapwith.siteno,  ]
		# remove swapwith.siteno from swapmat:
		swapmat <- swapmat[swapmat$site != swapwith.siteno,  ]
		n.swaps <- n.swaps + 1
		swap.listno <- NA
		if(any(site.val == 1))
			swap.listno <- min((1:n.sites)[site.val == 1])
	}
	# if there are any rows left in swapmat, the BAM procedure is not 
	# finished: need to form connected cliques and find chains, blocks, etc.

	list(swapmat = as.data.frame(swapmat), swap = n.swaps)
}
twoblock.func <- function(adjmat)
{
	# twoblock.func
	# Assume that all singletons and chains have been eliminated already.
	# The "s"-sites in the 2-block must therefore have valency exactly 2.
	#
	# adjmat is input with row-names and column-names equal to the 
	# SITE names of the corresponding sites.
	#
	# Search for 2-blocks with "r"-sites in the columns of adjmat: 
	# i.e. search for 2 or more rows (s-sites) with total valency contained 
	# within the same 2 columns (r-sites). 
	# (If no blocks are found, the adjmat matrix is transposed and sent
	# to this function again.)
	#
	s.valency <- apply(adjmat, 1, sum)
	workmat <- adjmat[s.valency == 2,  ]
	n.rows <- length(s.valency[s.valency==2])
	# workmat is the working matrix: only contains s-sites with valency exactly 2.
	# the "working numbers" of s-sites refer to their row position in workmat;
	# they are not connected to the real site numbers of these sites, which
	# are given by the rownames of workmat. 
	# Function returns the SITE NUMBERS of the block members. 
	if(n.rows<2) status <- "none" 
	else {		
		s.sites <- apply(workmat, 1, paste, collapse = " ")
		nwork.s <- length(s.sites)
		workblock.s <- (1:nwork.s)[duplicated(s.sites)]
		# nwork.s is the number of s-sites in the working matrix. 
		# workblock.s gives the working numbers of the s-sites that are
		# potentially in a 2-block.
		# 
		# if there are no 2-blocks, workblock.s will be numeric(0)
		# and have length=0.
		if(length(workblock.s) > 0) {
			select.s.site <- workblock.s[1]
			inblock.s <- s.sites[s.sites == s.sites[select.s.site]]
			# now reclaim the site numbers of the sites in the block, 
			# using the row-names and column-names of adjmat.
			inblock.s.sitenos <- names(inblock.s)
			select.s.siteno <- inblock.s.sitenos[1]
			inblock.r.sitenos <- dimnames(workmat)[[2]][workmat[
				select.s.siteno,  ] == 1]
			status <- "found"
			block <- list(r = as.numeric(inblock.r.sitenos), s = 
				as.numeric(inblock.s.sitenos))
		}
		else status <- "none"
	}
	if(status == "none")
		block <- NULL
	list(block = block, status = status)
}
rsblock.wrap <- function(mat, r)
{
	# rsblock.wrap: wrapper for rsblock.func
	rsblock.func(mat, r, r, NULL, NULL, as.numeric(dimnames(mat)[[2]]),
		as.numeric(dimnames(mat)[[1]]), F)
}
catline <- function(...)
{
	cat(..., "\n")
}