Use GBL v3 within hps-java

tracking/src/main/java/org/hps/recon/tracking/gbl/GBLexampleJna1.java

@@ -6,7 +6,14 @@
 import java.util.List;
 import java.util.ArrayList;
 import static java.lang.Math.sqrt;
+
 import org.apache.commons.math3.distribution.NormalDistribution;
+import org.apache.commons.cli.CommandLine;
+import org.apache.commons.cli.ParseException;
+import org.apache.commons.cli.HelpFormatter;
+import org.apache.commons.cli.DefaultParser;
+import org.apache.commons.cli.Options;
+
 import org.lcsim.util.aida.AIDA;
 
 import com.sun.jna.ptr.IntByReference;
@@ -15,16 +22,56 @@
 //import hep.aida.IHistogram1D;
 import java.io.IOException;
 
+/**
+ * java translation of GBL's example 1 using GBL via JNA
+ * <p>
+ * This translation is helpful for a few reasons.
+ * <ol>
+ *  <li>
+ *    During development, it can be used to check for basic
+ *    functionality without having to have a full input data
+ *    file.
+ *  </li>
+ *  <li>
+ *    Validation of developments can be done by using this
+ *    to compare two different branches without worrying about
+ *    whether the input data is affecting the comparison.
+ *  </li>
+ *  <li>
+ *    It is a helpful example for how the GBL-JNA classes
+ *    should be used in this java package.
+ *  </li>
+ * </ol>
+ * <p>
+ * <h2>Usage</h2>
+ * This class, while it is called in the unit testing,
+ * can also be called as its own executable. Since this
+ * example uses packages from outside the tracking package,
+ * it should be run from the central hps-distribution jar.
+ * <p>
+ * <code>
+ * java \
+ *   -cp distribution/target/hps-distribution-XXX-SNAPSHOT-bin.jar \
+ *   org.hps.recon.tracking.gbl.GBLexampleJna1
+ * </code>
+ */
 public class GBLexampleJna1 {
+
+    private NormalDistribution norm = new NormalDistribution();
 
     private int nTry = 100000;
     private int nLayer = 10;
-    private NormalDistribution norm = new NormalDistribution();
     private String outputPlots = "example1.root";
     private boolean debug = false;
 
-
     public AIDA aida;
+
+    public GBLexampleJna1(int nTry, int nLayer, boolean debug, String outputPlots) {
+        this.nTry = nTry;
+        this.nLayer = nLayer;
+        this.debug = debug;
+        this.outputPlots = outputPlots;
+    }
 
     public void setupPlots() {
         aida = AIDA.defaultInstance();
@@ -50,6 +97,7 @@ public void setupPlots() {
     public void runExample() {
         setupPlots();
         System.out.println("Running GBL Example!");
+        System.out.println("  N Tries = "+nTry+" with N Layers = "+nLayer);
         long startTime = System.nanoTime();
 
         double sinLambda = 0.3;
@@ -231,7 +279,7 @@ public void runExample() {
                     System.out.println("Points::" + listOfPoints.size());
                 }
 
-                jacPointToPoint = gblSimpleJacobianSvn(step, cosLambda, bfac);
+                jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac);
 
                 clPar = jacPointToPoint.times(clPar);
                 clCov = jacPointToPoint.times(clCov.times(jacPointToPoint.transpose()));
@@ -288,25 +336,19 @@ public void runExample() {
             Matrix seed = new SymMatrix(5);
             seed.set(0,0,100000);
 
-            //GblTrajectoryJna traj = new GblTrajectoryJna(listOfPoints,1,1,1);
-            GblTrajectoryJna traj = new GblTrajectoryJna(listOfPoints,1,seed,1,1,1);
+            //GblTrajectoryJna traj = new GblTrajectoryJna(listOfPoints,true,true,true);
+            GblTrajectoryJna traj = new GblTrajectoryJna(listOfPoints,1,seed,true,true,true);
 
-            if (traj.isValid() == 0 ) {
+            if (traj.isValid()) {
                 System.out.println("Example1: " + " Invalid GblTrajectory -> skip");
             }
-            double[] dVals = new double[2];
-            int [] iVals = new int[1];
-
-            DoubleByReference Chi2 = new DoubleByReference(0.);
-            DoubleByReference lostWeight = new DoubleByReference(0.);
-            IntByReference Ndf = new IntByReference(0);
-
-            //traj.fit(dVals,iVals,"");
-
-            traj.fit(Chi2,Ndf,lostWeight,"");
-            //Chi2Sum += dVals[0];
-            //NdfSum += iVals[0];
-            //LostSum += dVals[1];
+
+            DoubleByReference Chi2 = new DoubleByReference();
+            IntByReference Ndf = new IntByReference();
+            DoubleByReference lostWeight = new DoubleByReference();
+
+            traj.fit(Chi2, Ndf, lostWeight, "");
+
             Chi2Sum += Chi2.getValue();
             NdfSum += Ndf.getValue();
             LostSum += lostWeight.getValue();
@@ -315,6 +357,12 @@ public void runExample() {
             //aida.histogram1D("Chi2").fill(dVals[0]);
             //aida.histogram1D("Ndf").fill(iVals[0]);
             //aida.histogram1D("Chi2_Ndf").fill(dVals[0]/(double)iVals[0]);
+
+            // memory cleanup
+            traj.delete();
+            for (GblPointJna pt : listOfPoints) {
+                pt.delete();
+            }
         }
 
         long endTime = System.nanoTime();
@@ -324,7 +372,6 @@ public void runExample() {
         System.out.printf("Chi2/Ndf = %f \n", Chi2Sum / (double) NdfSum);
         System.out.printf("Tracks Fitted  %d \n", numFit);
 
-
         if (outputPlots != null) {
             try {
                 aida.saveAs(outputPlots);
@@ -335,9 +382,8 @@ public void runExample() {
         }
     }
 
-
-    private Matrix gblSimpleJacobianSvn(double ds, double cosl, double bfac) {
-
+    // manually translated to java from C++ from GBL CPP exampleUtil source
+    private Matrix gblSimpleJacobian(double ds, double cosl, double bfac) {
         Matrix jac = new Matrix(5,5);
         jac.UnitMatrix();
         jac.set(1, 0, -bfac * ds * cosl);
@@ -347,6 +393,43 @@ private Matrix gblSimpleJacobianSvn(double ds, double cosl, double bfac) {
         return jac;
 
     }   
+
+    public static void main(String[] args) {
+        Options options = new Options();
+        options.addOption("v", false, "verbose debugging printouts");
+        options.addOption("h", false, "print help message and exit");
+        options.addOption("o", true , "output file for histograms of residuals");
+        options.addOption("t", true , "number of tries to run (default: 100)");
+        options.addOption("l", true , "number of layers to track through (default: 10)");
+
+        boolean debug = false;
+        int nTries = 100;
+        int nLayers = 10;
+        String outputFile = "example1.root";
+
+        DefaultParser parser = new DefaultParser();
+        try {
+            CommandLine cli = parser.parse(options, args);
+            if (cli.hasOption("h")) {
+                HelpFormatter help = new HelpFormatter();
+                help.printHelp("GBLexampleJna1", options);
+                return;
+            }
+            debug = cli.hasOption("v");
+            outputFile = cli.getOptionValue("o", outputFile);
+            if (cli.hasOption("t")) {
+                nTries = Integer.parseInt(cli.getOptionValue("t"));
+            }
+            if (cli.hasOption("l")) {
+                nLayers = Integer.parseInt(cli.getOptionValue("l"));
+            }
+        } catch (ParseException exp) {
+            System.err.println(exp.getMessage());
+        }
+
+        GBLexampleJna1 eg = new GBLexampleJna1(nTries,nLayers,debug,outputFile);
+        eg.runExample();
+    }
 }
 
 

tracking/src/main/java/org/hps/recon/tracking/gbl/GblHelixPrediction.java

@@ -1,34 +1,18 @@
 package org.hps.recon.tracking.gbl;
 
-import com.sun.jna.Library;
-import com.sun.jna.Native;
-import com.sun.jna.Pointer;
+import com.sun.jna.Pointer; 
 
 import org.hps.recon.tracking.gbl.matrix.Matrix;
 import org.hps.recon.tracking.gbl.matrix.Vector;
 
+/**
+ * wrapper class for GblHelixPrediction JNA functions
+ * <p>
+ * Re-promote the JNA functions into class member functions while
+ * adding some helpful translation.
+ */
 public class GblHelixPrediction {
 
-    public interface GblHelixPredictionInterface extends Library {
-
-        GblHelixPredictionInterface INSTANCE = (GblHelixPredictionInterface) Native.loadLibrary("GBL", GblHelixPredictionInterface.class);
-
-        Pointer GblHelixPredictionCtor(double sArc, double[] aPred,
-                                       double [] tDir, double [] uDir, double [] vDir,
-                                       double [] nDir, double [] aPos);
-
-        double GblHelixPrediction_getArcLength(Pointer self);
-        void GblHelixPrediction_getMeasPred(Pointer self, double [] prediction);
-        void GblHelixPrediction_getPosition(Pointer self, double[] position);
-        void GblHelixPrediction_getDirection(Pointer self,double[] direction);
-
-        double GblHelixPrediction_getCosIncidence(Pointer self);
-
-        //array for the curvilinear directions (2x3 matrix)
-        void GblHelixPrediction_getCurvilinearDirs(Pointer self, double [] curvilinear);
-
-    }
-
     private Pointer self;
 
     //aPred is a 2-vector
@@ -41,7 +25,7 @@ public GblHelixPrediction(double sArc, double[] aPred,
                               double [] tDir, double [] uDir, double [] vDir, 
                               double [] nDir, double [] aPos) {
 
-        self = GblHelixPredictionInterface.INSTANCE.GblHelixPredictionCtor(sArc, aPred,
+        self = GblInterface.INSTANCE.GblHelixPredictionCtor(sArc, aPred,
                                                                            tDir, uDir, vDir,
                                                                            nDir,aPos);
     }
@@ -50,6 +34,10 @@ public GblHelixPrediction(Pointer ptr) {
         self = ptr;
     }
 
+    public void delete() {
+        GblInterface.INSTANCE.GblHelixPrediction_delete(self);
+    }
+
     public GblHelixPrediction() {
 
         double sArc = 0;
@@ -60,7 +48,7 @@ public GblHelixPrediction() {
         double [] nDir = new double[2];
         double [] aPos = new double[2];
 
-        self = GblHelixPredictionInterface.INSTANCE.GblHelixPredictionCtor(sArc, aPred,
+        self = GblInterface.INSTANCE.GblHelixPredictionCtor(sArc, aPred,
                                                                            tDir, uDir, vDir,
                                                                            nDir,aPos);
 
@@ -76,37 +64,37 @@ public GblHelixPrediction(double sArc, Vector v_aPred, Vector v_tDir, Vector v_u
         double[] aPos  =  v_aPos.getColumnPackedCopy();
 
 
-        self = GblHelixPredictionInterface.INSTANCE.GblHelixPredictionCtor(sArc, aPred,
+        self = GblInterface.INSTANCE.GblHelixPredictionCtor(sArc, aPred,
                                                                            tDir, uDir, vDir,
                                                                            nDir, aPos);
     }
 
     public double getArcLength() {
-        return GblHelixPredictionInterface.INSTANCE.GblHelixPrediction_getArcLength(self);
+        return GblInterface.INSTANCE.GblHelixPrediction_getArcLength(self);
     }
 
     public void getMeasPred(double [] prediction) {
-        GblHelixPredictionInterface.INSTANCE.GblHelixPrediction_getMeasPred(self,prediction);
+        GblInterface.INSTANCE.GblHelixPrediction_getMeasPred(self,prediction);
     }
 
     public void getPosition(double [] position) {
-        GblHelixPredictionInterface.INSTANCE.GblHelixPrediction_getPosition(self, position);
+        GblInterface.INSTANCE.GblHelixPrediction_getPosition(self, position);
     }
 
     public void getDirection(double [] direction) {
-        GblHelixPredictionInterface.INSTANCE.GblHelixPrediction_getDirection(self, direction);
+        GblInterface.INSTANCE.GblHelixPrediction_getDirection(self, direction);
     }
 
     public double getCosIncidence() {
-        return GblHelixPredictionInterface.INSTANCE.GblHelixPrediction_getCosIncidence(self);
+        return GblInterface.INSTANCE.GblHelixPrediction_getCosIncidence(self);
     }
 
     public Matrix getCurvilinearDirs() {
 
         Matrix curvilinearDirs = new Matrix(2,3);
         double[] dirArray = new double[6];
 
-        GblHelixPredictionInterface.INSTANCE.GblHelixPrediction_getCurvilinearDirs(self, dirArray);
+        GblInterface.INSTANCE.GblHelixPrediction_getCurvilinearDirs(self, dirArray);
 
         curvilinearDirs.set(0,0,dirArray[0]);
         curvilinearDirs.set(0,1,dirArray[1]);