Weiszfeld method for the multi-facility location problem

MultiWeiszfeld.scala

/** 
 * 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)
  }
}