Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Weiszfeld method for the multi-facility location problem
/**
* A straightforward implementation of the multi-facility Weiszfeld method described by
* Iyigun et al. for 2D euclidean. http://dx.doi.org/10.1016/j.orl.2009.11.005
* Have some failings, e.g. if all the initial locations are the same
*
* @param points The demand points
* @param initial Initial locations of facilities
*/
def multiWeiszfeld(points: Seq[Point], initial: Seq[Point]): Seq[Point] = {
val eps = 1e-6
val K = initial.length
val px = points.map(_.x).toArray
val py = points.map(_.y).toArray
var diff = 0.0
var z0, z1 = initial
do {
val z1x = z1.map(_.x).toArray
val z1y = z1.map(_.y).toArray
val newCenters_x = (0 until K).map(k => recenter(px, z1x, k))
val newCenters_y = (0 until K).map(k => recenter(py, z1y, k))
val newCenters = newCenters_x.zip(newCenters_y).map(t => Point(t._1, t._2))
z0 = z1
z1 = newCenters
diff = (0 until z0.size).map(i => z0(i).distance(z1(i))).sum
} while (z0 != z1 && diff > eps)
z1
}
/**
* @param points The demand points
* @param K number of facilities
*/
def multiWeiszfeld(points: Seq[Point], K: Int): Seq[Point] = {
multiWeiszfeld(points, points.take(K))
}
def prob(x: Array[Double], c: Array[Double], k: Int, i: Int) = {
val K = c.length
val num = ((0 until k) ++ (k + 1 until K)).map { j =>
math.abs(x(i) - c(j))
}.reduceLeft(_ * _)
val den = (0 until K).map { l =>
((0 until l) ++ (l + 1 until K)).map { m =>
math.abs(x(i) - c(m))
}.reduceLeft(_ * _)
}.sum
if (num / den isNaN) 0 else num / den
}
def recenter(x: Array[Double], c: Array[Double], k: Int) = {
val N = x.length
val w = 1.0
val ret = (0 until N).map { i =>
val num = (prob(x, c, k, i) * w) / math.abs(x(i) - c(k))
val den = (0 until N).map { j =>
val v = (prob(x, c, k, j) * w) / math.abs(x(j) - c(k))
if (v.isNaN() || v.isInfinite()) 0 else v
}.sum
if (num.isInfinite || num.isNaN) 0 else (num / den) * x(i)
}.sum
ret
}
case class Point(x: Double, y: Double) {
def distance(that: Point) = {
val x = this.x - that.x
val y = this.y - that.y
math.sqrt(x*x + y*y)
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.