Daily Programmer 321 (Reashu)

import org.apache.commons.math.linear.Array2DRowRealMatrix
import org.apache.commons.math.linear.ArrayRealVector
import org.apache.commons.math.linear.LUDecompositionImpl
import org.apache.commons.math.linear.RealVector

/**
 * https://www.reddit.com/r/dailyprogrammer/comments/6ksmh5/20170630_challenge_321_hard_circle_splitter/
 */
val EPSILON = 0.00001

fun main(args: Array<String>) {
    val inputCoordinates = args.drop(1).map(String::toDouble)
    val xs = inputCoordinates.filterIndexed { index, _ -> even(index) }
    val ys = inputCoordinates.filterIndexed { index, _ -> odd(index) }
    val pairedCoordinates = xs.zip(ys);
    val points = pairedCoordinates.map { PointX(arrayOf(it.first, it.second)) }

    val bestCircle = solve(points)

    // TODO Iterate on the circle, increasing / decreasing size slightly to compensate for the errors in matrix solving

    if (bestCircle != null) {
        println(bestCircle.center.vector.data.map(Double::toString).joinToString(separator = " "))
        println(bestCircle.radius)
    } else {
        println("No solution")
    }
}

fun even(i: Int) = i % 2 == 0
fun odd(i: Int) = i % 2 == 1

fun solve(points: List<PointX>): CircleX? {
    val scoredPoints = points.map({ point ->
        val closestPoints = closestHalf(point, points)
        val score = 1 / closestPoints.last().euclideanDistanceTo(point)
        Triple(point, score, closestPoints)
    }).sortedByDescending{ it.second }

    // Create an n-dimensional "square" at 0.5, 0.5, ...
    val squareToFit = SquareX(PointX(Array(points[0].dimensions, {0.5})), 1.0)
    for (point in scoredPoints) {
        val minimumCircle = minimumBallMTF(point.third, listOf())
        if (squareToFit.contains(minimumCircle)) {
            return minimumCircle
        }
    }
    return null
}

/**
 * Finds the N/2 points closest to this point among the given N points, sorted by ascending distance
 */
fun closestHalf(primary: PointX, all: List<PointX>): List<PointX> {
    val sortedPoints = all.sortedBy { it.euclideanDistanceTo(primary) }
    return sortedPoints.subList(0, all.size / 2)
}

data class PointX(val vector: RealVector) {
    val dimensions = vector.dimension

    constructor(positions: Array<Double>): this(ArrayRealVector(positions))

    fun euclideanDistanceTo(other: PointX): Double = vector.getDistance(other.vector)
    fun longestOrthogonalDistanceTo(other: PointX): Double = vector.getLInfDistance(other.vector)
}

fun origo(dimensions: Int): PointX {
    return PointX(Array(dimensions, {0.0}))
}

data class CircleX(val center: PointX, val radius: Double) {
    fun contains(point: PointX): Boolean {
        // Sometimes we need an "empty" circle without knowing how many dimensions it has
        // That circle should never contain anything
        if (center.dimensions == 0) return false

        return center.euclideanDistanceTo(point) <= radius + EPSILON;
    }
    val area = Math.pow(radius, center.dimensions.toDouble())
}

data class SquareX(val center: PointX, val sideLength: Double) {
    fun contains(point: PointX): Boolean {
        return center.longestOrthogonalDistanceTo(point) <= sideLength / 2 + EPSILON
    }

    fun contains(circle: CircleX): Boolean {
        return center.longestOrthogonalDistanceTo(circle.center) + circle.radius <= sideLength / 2 + EPSILON
    }
}


// Based on https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf
fun minimumBallMTF(points: List<PointX>, boundaryPoints: List<PointX>): CircleX {
    var l = points;
    var b = boundaryPoints;

    // Find the smallest ball which fits our boundary points
    var mb = ballForBoundary(b)
    if (l.isEmpty()) return mb

    // If we have enough boundary points, it is guaranteed to fit all points
    if (b.size == points[0].dimensions + 1)
        return mb
    // If not, find the points which do not fit and calculate a new ball which includes them
    for (i in 0 until l.size) {
        if (!mb.contains(l[i])) {
            // Calculate a new ball with all previously included points, and the failing point on the edge
            mb = minimumBallMTF(l.subList(0, i), b.plus(l[i]))
            /*
             * Reorder the list of points so that the now-included ones are
             *  (1) kept track of and
             *  (2) checked early in future recursive calls, as they have proved "important"
             * Reordering the beginning of the list is safe because all of those points are included, and
             * iteration will continue with the next potentially unincluded point
             */
            l = listOf(l[i]).union(l).toList()
        }
    }
    return mb
}

fun ballForBoundary(boundaryPoints: List<PointX>): CircleX {
    if (boundaryPoints.size == 0)
        return CircleX(PointX(arrayOf()), 0.0)
    if (boundaryPoints.size == 1)
        return CircleX(boundaryPoints[0], 0.0)

    val firstBoundaryVector = boundaryPoints[0].vector
    val relativeBoundaryVectors = boundaryPoints.drop(1).map({ it.vector.subtract(firstBoundaryVector) })

    val a = Array2DRowRealMatrix(relativeBoundaryVectors.size, relativeBoundaryVectors.size)
    for (row in 0 until relativeBoundaryVectors.size) {
        for (col in 0 until relativeBoundaryVectors.size) {
            a.setEntry(row, col, 2.0 * (relativeBoundaryVectors[row].dotProduct(relativeBoundaryVectors[col])))
        }
    }
    val b = ArrayRealVector(relativeBoundaryVectors.map({ it.dotProduct(it) }).toDoubleArray())
    val solver = LUDecompositionImpl(a).solver
    try {
        val x = solver.solve(b)
        val dimensions = boundaryPoints[0].dimensions
        val relativeCenter = x.data.zip(relativeBoundaryVectors).fold(origo(dimensions).vector, { acc, (weight, point) ->
            acc.add(point.mapMultiply(weight))
        })
        val radius = Math.sqrt(relativeCenter.dotProduct(relativeCenter))
        return CircleX(PointX(relativeCenter.add(firstBoundaryVector)), radius);
    } catch (e: Exception) {
        throw RuntimeException("Error when solving for points $boundaryPoints: ${e.message}")
    }
}