Skip to content

Commit

Permalink
Merge pull request #958 from tomeichlersmith/gbl-v3-jna
Browse files Browse the repository at this point in the history
Use GBL v3 within hps-java
  • Loading branch information
cbravo135 authored Apr 20, 2023
2 parents e0695c4 + 4f599a3 commit 585058f
Showing 13 changed files with 515 additions and 318 deletions.
129 changes: 106 additions & 23 deletions tracking/src/main/java/org/hps/recon/tracking/gbl/GBLexampleJna1.java
Original file line number Diff line number Diff line change
@@ -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();
}
}


Original file line number Diff line number Diff line change
@@ -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]);
Loading

0 comments on commit 585058f

Please sign in to comment.