STEPS, Institute of Transportation Studies, Davis
Personal tools
Sections
Navigation
STEPS, Institute of Transportation Studies, Davis

2028 Academic Surge
University of California, Davis
One Shields Avenue
Davis, CA 95616

(530) 752-6548 phone

 

Weiszfeld's Algorithm in R

R code to run Weiszfeld's Algorithm. This algorithm finds a central tendency point for a set of weighted points. It includes a simple animation of the algorithms movement towards the central point.

#Weiszfeld's Algorithm in [R]

# See Francis, R.L., et al, "Facility Layout and Location: an analytical approach",
# Prentice-Hall for algorithm reference

# Find the center of mass for a network of weighted points

# INPUT: a data frame (a,b,w) = a (x) coordinate, b (y) coordinate, weight
#		 iters = number of iterations (default to 25)
#        start = search start point (will default)
#   Example weiszfeld(x,iters=20,start=c(a,b))

# OUTPUT: matrix of x,y,f estimates
# The f is the minimized sum, x and y are estimate coordinates. 

weiszfeld = function(dat=NULL, iters=25, start=NULL,speed=0.1){
	
	if(is.null(dat)){
		### REMOVE this block
		cat("Null data. Running in example mode.\n")
		mat = matrix(c(2,5,4,1,3,7,8,8,1,1,1,2),nrow=4,ncol=3)
		dat = as.data.frame(mat)
		names(dat)=c("a","b","w")
		
	}
	

	# check dat col names
	cols = c("a","b","w")
	for(i in 1:3){
		if(cols[i] != names(dat)[i]){
			stop("Column names for input dataframe must be: a (x coord), b (y coord), w (weight)")
		}
	}
	if(is.null(start)){
		# guess a start value
		start=c(mean(dat$a),mean(dat$b))
	}
	# Make sure you're not starting on top an existing point
	go=TRUE
	while(go==TRUE){
		# Make sure the start point does not duplicate a point in the list.
		o = subset(dat,a==start[1] && b==start[2])
		if(nrow(o)==0){
			go=FALSE
		} else {
			start[1]=start[1]+0.01
			start[2]=start[2]+0.01
		}
	}

	gamma.i = function(x,y,a,b,w){
		out = w/sqrt((x-a)^2+(y-b)^2)
		return(out)
	}
	
	cap.gamma = function(dat,start ){
		out = 0
		for(i in 1:nrow(dat)){
			row = dat[i,]
			out = out + gamma.i(start[1],start[2],row$a,row$b,row$w)
		}
		return(out)
			
	}
	
	calc.f = function(start,dat){
		out = 0
		for(i in 1:nrow(dat)){
			row = dat[i,]
			out = out + (row$w*sqrt((start[1]-row$a)^2+(start[2]-row$b)^2))
		}
		return(out)
	}
	
	ret=data.frame(matrix(0,nrow=iters,ncol=3))
	names(ret)=c("x","y","f")
	j = 1
	while(iters > 0){
		print(start)
		plot(dat$a,dat$b)
		points(start[1],start[2],pch="x")
		Gamma = cap.gamma(dat,start)
		new.x = 0
		new.y = 0
		for (i in 1:nrow(dat)){
			row = dat[i,]
			new.x=new.x + ((gamma.i(start[1],start[2],row$a,row$b,row$w)*row$a)/Gamma)
			new.y=new.y + ((gamma.i(start[1],start[2],row$a,row$b,row$w)*row$b)/Gamma)

		}
		# calculate f
		
		ret[j,]=c(start[1], start[2], (calc.f(start,dat)))
		start[1]=new.x
		start[2]=new.y
		iters=iters-1
		j = j+1
		Sys.sleep(speed)
	}
	return(ret)
}