
/* UPMC - c.durr - 2016-2017
   Conception et Pratique de l'Algorithmique
   TME 4 - algorithme k-means++
*/

import java.util.*;
import java.io.*;

class Point {
    double lat, lon;

    Point(double lat, double lon) {
        this.lat = lat;
        this.lon = lon;
    }

    double dist(Point p) {
        double dlon = p.lon - lon;
        double dlat = p.lat - lat;
        double a = Math.pow(Math.sin(Math.toRadians(dlat / 2)), 2) +
                   Math.cos(Math.toRadians(lat)) * Math.cos(Math.toRadians(p.lat)) *
                    Math.pow(Math.sin(Math.toRadians(dlon / 2)), 2);
        return Math.asin( Math.sqrt(a) ) * 6373;
    }
}

class Kmeansplusplus {

    static ArrayList<Point> lire(String filename) throws FileNotFoundException {
        // lit le fichier d'entrée, le recopie dans la sortie et retourne la liste de ses points
        // complexité: linéaire
        ArrayList<Point> points = new ArrayList<Point>();
        Scanner in = new Scanner(new File(filename));
        in.useLocale(Locale.US); // pour utiliser le point comme séparateur des décimaux, et pas la virgule
        while (in.hasNext()) {
            String v = in.next();  // consomme le 'v'
            assert v == "v";
            long id = in.nextLong();
            double lat = in.nextDouble();
            double lon = in.nextDouble();
            //                   -- recopier sur la sortie
            System.out.format(Locale.US, "v %d %.10f %.10f\n",id, lat, lon);
            points.add(new Point(lat, lon));
        }
        return points;
    }

    static Point [] init(ArrayList<Point> points, int k) {
        // produit une liste de k centres choisit selon l'algorithme k-means++
        // complexité: O(nk^2) pour n=points.length()
        Random alea = new Random();
        Point centers[] = new Point[k];
        centers[0] = points.get(alea.nextInt(points.size()));
        for (int j=1; j < k; j++) {
            double total = 0.0;
            for (Point p: points) {
                double d = Double.MAX_VALUE;
                for (int i=0; i < j; i++) {
                    d = Math.min(d, p.dist(centers[i]));
                }
                double d2 = d*d;
                total += d2;
                if (alea.nextDouble() * total < d2)   // ne pas diviser par total, pourrait être 0
                    centers[j] = p;
            }
        }
        return centers;
    }

    static double Lloyd(ArrayList<Point> points, Point centers[]) {
        // amélioration itérative de Lloyd, retourne la valeur objective finale atteinte au point fixe, modifie centers
        // complexité: O(nkI) où I est le nombre d'itérations
        int k = centers.length;
        int n = points.size();
        double previous_obj = Double.MAX_VALUE, obj =  Double.MAX_VALUE / 2;
        while (previous_obj > obj) {
            previous_obj = obj;
            obj = 0.0;                                  // recalculer l'objectif
            double tot_lat [] = new double[k];          // calculer les centres de gravité
            double tot_lon [] = new double[k];
            int tot_siz [] = new int[k];                // taille de chaque cluster
            for (Point p : points) {                    // trouver le centre le plus proche de p
                int closest_center = -1;                // rien pour l'instant
                double dist_center = Double.MAX_VALUE;
                for (int j=0; j < k; j++) {
                    double d = p.dist(centers[j]);
                    if (d < dist_center) {
                        dist_center = d;
                        closest_center = j;
                    }
                }
                tot_lat[closest_center] += p.lat;       // cumuler les coordonnées pour déterminer la moyenne
                tot_lon[closest_center] += p.lon;
                tot_siz[closest_center] += 1;
                obj += dist_center * dist_center;
            }
            for (int j=0; j < k; j++)
                if (tot_siz[j] > 0) {
                    centers[j].lat = tot_lat[j] / tot_siz[j];  // affecter aux centres de gravité
                    centers[j].lon = tot_lon[j] / tot_siz[j];
                }
            System.err.println("objective value = " + obj);    // montrer progression
        }
        return obj;
    }

    public static void main(String[] args) throws FileNotFoundException {
        int k = Integer.parseInt(args[0]);
        ArrayList<Point> points = lire(args[1]);        // lire l'entrée
        Point centers [] = init(points, k);            // choisir k centers avec l'algo k-means++
        double obj = Lloyd(points, centers);            // recherche locale par l'algorithme de Lloyd
        for (Point p: centers)                          // afficher la solution
            System.out.format(Locale.US, "s %.10f %.10f\n", p.lat, p.lon);
        System.out.format(Locale.US, "o %.10g\n", obj);
    }
}
