Fixes to logP and cov functions for GMMs. (#133)

Co-authored-by: samypr100 <3933065+samypr100@users.noreply.github.com>

Author: Birdasaur

build.gradle

@@ -367,7 +367,7 @@ assemble {
 jlink {
     def fileSep = System.getProperty('file.separator')
     def imageZipFile = layout.buildDirectory.file("${artifactNameUpper}-${project.version}.zip")
-    options.set(['--strip-debug', '--compress', '2', '--no-header-files', '--no-man-pages'])
+    options.set(['--strip-debug', '--compress', 'zip-6', '--no-header-files', '--no-man-pages'])
     imageZip.set(imageZipFile)
     launcher {
         def currentOS = org.gradle.internal.os.OperatingSystem.current()

src/main/java/edu/jhuapl/trinity/javafx/components/panes/Shape3DControlPane.java

@@ -121,7 +121,7 @@ private void buildFindClustersTab() {
         findClustersTab.setContent(findClusterBorderPane);
 
         componentsSpinner = new Spinner(
-            new SpinnerValueFactory.IntegerSpinnerValueFactory(2, 20, 5, 1));
+            new SpinnerValueFactory.IntegerSpinnerValueFactory(2, 500, 5, 1));
         componentsSpinner.setPrefWidth(SPINNER_PREF_WIDTH);
         componentsSpinner.setEditable(true);
         iterationsSpinner = new Spinner(

src/main/java/edu/jhuapl/trinity/javafx/handlers/ManifoldEventHandler.java

@@ -298,11 +298,6 @@ public void handleExport(ManifoldEvent event) {
     }
 
     public void handleNewManifoldData(ManifoldEvent event) {
-//        Platform.runLater(() -> {
-//            App.getAppScene().getRoot().fireEvent(
-//                new CommandTerminalEvent("Loading Manifold Data...",
-//                    new Font("Consolas", 20), Color.GREEN));
-//        });
 //        System.out.println("Loading Manifold Data...");
         ManifoldData md = (ManifoldData) event.object1;
         //convert deserialized points to Fxyz3D point3ds

src/main/java/edu/jhuapl/trinity/javafx/javafx3d/AsteroidFieldPane.java

@@ -757,7 +757,7 @@ private void resetAsteroids() {
 
     public Manifold3D makeHull(List<Point3D> labelMatchedPoints, String label, Double tolerance) {
         Manifold3D manifold3D = new Manifold3D(
-            labelMatchedPoints, true, true, true, tolerance
+            labelMatchedPoints, true, false, false, tolerance
         );
         manifold3D.quickhullMeshView.setCullFace(CullFace.FRONT);
 //        manifold3D.addEventHandler(MouseEvent.MOUSE_CLICKED, e -> {

src/main/java/edu/jhuapl/trinity/javafx/javafx3d/Manifold3D.java

@@ -77,7 +77,7 @@ public class Manifold3D extends Group {
 
     public Manifold3D(List<Point3D> point3DList, boolean triangulate, boolean makeLines, boolean makePoints, Double tolerance) {
         originalPoint3Ds = point3DList;
-        buildHullMesh(point3DList, triangulate, makeLines, makePoints, tolerance);
+        buildHullMesh(point3DList, triangulate, tolerance);
 
         List<Point3D> fxyzPoints = new ArrayList<>();
         for (int i = 0; i < hull.getNumVertices(); i++) {
@@ -349,7 +349,7 @@ public void refreshMesh(List<Point3D> point3DList, boolean triangulate, boolean
         quickhullLinesTriangleMesh.getPoints().clear();
         quickhullLinesTriangleMesh.getTexCoords().clear();
         quickhullLinesTriangleMesh.getFaces().clear();
-        buildHullMesh(point3DList, triangulate, makeLines, makePoints, tolerance);
+        buildHullMesh(point3DList, triangulate, tolerance);
         quickhullMeshView.setMesh(quickhullTriangleMesh);
         if (makeLines) {
             quickhullLinesTriangleMesh.getPoints().addAll(quickhullTriangleMesh.getPoints());
@@ -361,7 +361,7 @@ public void refreshMesh(List<Point3D> point3DList, boolean triangulate, boolean
 //            makeDebugPoints(hull, artScale, false);
     }
 
-    private void buildHullMesh(List<Point3D> point3DList, boolean triangulate, boolean makeLines, boolean makePoints, Double tolerance) {
+    private void buildHullMesh(List<Point3D> point3DList, boolean triangulate, Double tolerance) {
         hull = new QuickHull3D();
         if (null != tolerance)
             hull.setExplicitDistanceTolerance(tolerance);
@@ -455,19 +455,21 @@ public void handle(long now) {
     }
 
     public void makeLines() {
+        boolean wasVisible = null != quickhullLinesMeshView
+            ? quickhullLinesMeshView.isVisible() : false;
         quickhullLinesTriangleMesh = new TriangleMesh();
         quickhullLinesTriangleMesh.getPoints().addAll(quickhullTriangleMesh.getPoints());
         quickhullLinesTriangleMesh.getTexCoords().addAll(quickhullTriangleMesh.getTexCoords());
         quickhullLinesTriangleMesh.getFaces().addAll(quickhullTriangleMesh.getFaces());
 
         quickhullLinesMeshView = new MeshView(quickhullLinesTriangleMesh);
-        PhongMaterial quickhullLinesMaterial = new PhongMaterial(Color.BLUE);
-        quickhullLinesMaterial.setSpecularColor(Color.BLUE); //fix for aarch64 Mac Ventura
+        PhongMaterial quickhullLinesMaterial = new PhongMaterial(Color.ALICEBLUE);
+        quickhullLinesMaterial.setSpecularColor(Color.ALICEBLUE); //fix for aarch64 Mac Ventura
         quickhullLinesMeshView.setMaterial(quickhullLinesMaterial);
         quickhullLinesMeshView.setDrawMode(DrawMode.LINE);
         quickhullLinesMeshView.setCullFace(CullFace.NONE);
         quickhullLinesMeshView.setMouseTransparent(true);
-
+        quickhullLinesMeshView.setVisible(wasVisible);
         getChildren().add(quickhullLinesMeshView);
     }
 
@@ -489,14 +491,15 @@ public void makeDebugPoints(QuickHull3D hull, float scale, boolean print) {
                     sb.append(", ");
             }
 
-            Sphere sphere = new Sphere(2.5);
-            PhongMaterial mat = new PhongMaterial(Color.BLUE);
-            mat.setSpecularColor(Color.BLUE); // fix for aarch64 Mac Ventura
+            Sphere sphere = new Sphere(1.5);
+            PhongMaterial mat = new PhongMaterial(Color.ALICEBLUE);
+            mat.setSpecularColor(Color.ALICEBLUE); // fix for aarch64 Mac Ventura
             sphere.setMaterial(mat);
             sphere.setTranslateX(point3D.x);
             sphere.setTranslateY(point3D.y);
             sphere.setTranslateZ(point3D.z);
             extrasGroup.getChildren().add(sphere);
+            sphere.setVisible(false);
 
             Label newLabel = new Label(String.valueOf(i));
             labelGroup.getChildren().addAll(newLabel);

src/main/java/edu/jhuapl/trinity/utils/clustering/ClusterUtils.java

@@ -3,6 +3,9 @@
 import org.slf4j.Logger;
 import org.slf4j.LoggerFactory;
 
+import java.util.ArrayList;
+import java.util.List;
+
 import static java.lang.Math.abs;
 import static java.lang.Math.sqrt;
 
@@ -56,6 +59,29 @@ public static double squaredDistanceWithMissingValues(double[] x, double[] y) {
         return dist;
     }
 
+    public static List<List<double[]>> extractGMMClusters(
+        double[][] data,
+        GaussianMixture gmm,
+        double threshold) {
+        List<List<double[]>> clusterPoints = new ArrayList<>();
+        for (int i = 0; i < gmm.components.length; i++) {
+            clusterPoints.add(new ArrayList<>());
+        }
+
+        for (double[] x : data) {
+            double[] post = gmm.posteriori(x);
+            int k = ClusterUtils.whichMax(post);
+            if (post[k] >= threshold) {
+                clusterPoints.get(k).add(x);
+            }
+        }
+
+        // Filter out tiny/degenerate clusters if needed
+        clusterPoints.removeIf(list -> list.size() < 4);
+
+        return clusterPoints;
+    }
+
     /**
      * Returns the sum of an array.
      *

src/main/java/edu/jhuapl/trinity/utils/clustering/GaussianDistribution.java

@@ -200,18 +200,20 @@ public double scatter() {
         return sigmaDet;
     }
 
-    public double logp(double[] x) {
-        if (x.length != dim) {
-            throw new IllegalArgumentException("Sample has different dimension.");
-        }
-
+    public double mahalanobis2(double[] x) {
         double[] v = x.clone();
         ClusterUtils.sub(v, mu);
-//        double result = sigmaInv.xAx(v) / -2.0;
-//        double[] Ax = mv(x);
-        double[] Ax = sigmaInv.operate(v);
-        double result = ClusterUtils.dot(x, Ax) / -2.0;
-        return result - pdfConstant;
+        double[] Av = sigmaInv.operate(v);
+        return ClusterUtils.dot(v, Av);
+    }
+
+    public double logp(double[] x) {
+        if (x.length != dim) throw new IllegalArgumentException("Sample has different dimension.");
+        double[] v = x.clone();
+        ClusterUtils.sub(v, mu);                 // v = x - μ
+        double[] Av = sigmaInv.operate(v);       // Σ⁻¹ v
+        double quad = ClusterUtils.dot(v, Av);   // vᵀ Σ⁻¹ v
+        return -0.5 * quad - pdfConstant;
     }
 
     public double p(double[] x) {
@@ -472,6 +474,10 @@ public double logLikelihood(double[][] x) {
         return L;
     }
 
+    public int dim() {
+        return dim;
+    }
+
     @Override
     public String toString() {
         return String.format("Gaussian(mu = %s, sigma = %s)", Arrays.toString(mu), sigma);

src/main/java/edu/jhuapl/trinity/utils/clustering/GaussianMixture.java

@@ -240,30 +240,45 @@ public double[] mean() {
     }
 
     public RealMatrix cov() {
-        double w = components[0].priori();
-        RealMatrix v = components[0].distribution().cov();
-
-        int m = v.getRowDimension();
-        int n = v.getColumnDimension();
-        RealMatrix cov = MatrixUtils.createRealMatrix(m, n);
+        double[] mu = mean();
+        int d = mu.length;
+        RealMatrix C = MatrixUtils.createRealMatrix(d, d);
 
-        for (int i = 0; i < m; i++) {
-            for (int j = 0; j < n; j++) {
-                cov.setEntry(i, j, w * w * v.getEntry(i, j));
-            }
+        // within-component variance
+        for (GaussianMixtureComponent c : components) {
+            double w = c.priori();
+            RealMatrix Sk = c.distribution().cov();
+            C = C.add(Sk.scalarMultiply(w));
         }
 
-        for (int k = 1; k < components.length; k++) {
-            w = components[k].priori();
-            v = components[k].distribution().cov();
-            for (int i = 0; i < m; i++) {
-                for (int j = 0; j < n; j++) {
-                    cov.addToEntry(i, j, w * w * v.getEntry(i, j));
-                }
-            }
+        // between-component variance
+        for (GaussianMixtureComponent c : components) {
+            double w = c.priori();
+            double[] mk = c.distribution().mean();
+            double[] diff = mk.clone();
+            ClusterUtils.sub(diff, mu);
+            RealMatrix outer = MatrixUtils.createColumnRealMatrix(diff)
+                .multiply(MatrixUtils.createRowRealMatrix(diff));
+            C = C.add(outer.scalarMultiply(w));
         }
+        return C;
+    }
+
+    public boolean inDistribution(double[] x, double q) {
+        // choose most responsible component
+        double[] r = posteriori(x);
+        int idx = ClusterUtils.whichMax(r);
+        GaussianMixtureComponent c = components[idx];
+        double d2 = c.distribution().mahalanobis2(x);
+        // chi-square threshold
+        org.apache.commons.math3.distribution.ChiSquaredDistribution chi =
+            new org.apache.commons.math3.distribution.ChiSquaredDistribution(c.distribution().dim());
+        double thresh = chi.inverseCumulativeProbability(q);
+        return d2 <= thresh;
+    }
 
-        return cov;
+    public boolean inDistributionByLogP(double[] x, double tau) {
+        return Math.log(p(x)) >= tau;
     }
 
     public Pair<Integer, Double> maxPostProb(double[] x) {

src/main/resources/edu/jhuapl/trinity/fxml/ManifoldControl.fxml

@@ -450,9 +450,9 @@
                                                         </HBox>
                                                         <HBox alignment="CENTER" spacing="10.0" GridPane.rowIndex="2">
                                                             <children>
-                                                                <CheckBox fx:id="showWireframeCheckBox" mnemonicParsing="false" selected="true"
+                                                                <CheckBox fx:id="showWireframeCheckBox" mnemonicParsing="false" selected="false"
                                                                           text="Show Wire Frame"/>
-                                                                <CheckBox fx:id="showControlPointsCheckBox" mnemonicParsing="false" selected="true"
+                                                                <CheckBox fx:id="showControlPointsCheckBox" mnemonicParsing="false" selected="false"
                                                                           text="Show Control Points"/>
                                                             </children>
                                                         </HBox>