Commit in hps-java on MAIN | |||
pom.xml | +6 | -1 | 1.95 -> 1.96 |
src/test/java/org/lcsim/hps/recon/tracking/TestRunTrackReconTest.java | +47 | -54 | 1.5 -> 1.6 |
src/main/java/org/lcsim/hps/util/CompareHistograms.java | +40 | -22 | 1.1 -> 1.2 |
+93 | -77 |
add kolmogorov-smirnov test
diff -u -r1.95 -r1.96 --- pom.xml 4 Jun 2013 21:07:59 -0000 1.95 +++ pom.xml 3 Aug 2013 00:57:03 -0000 1.96 @@ -1,6 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?> <project xmlns="http://maven.apache.org/POM/4.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
- xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">
+ xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">
<modelVersion>4.0.0</modelVersion> <artifactId>hps-java</artifactId> <groupId>org.lcsim</groupId>
@@ -275,5 +275,10 @@
<artifactId>freehep-jaida-remote</artifactId> <version>3.3.0-12</version> </dependency>
+ <dependency> + <groupId>org.apache.commons</groupId> + <artifactId>commons-math3</artifactId> + <version>3.2</version> + </dependency>
</dependencies> </project>
diff -u -r1.5 -r1.6 --- TestRunTrackReconTest.java 13 Jul 2013 01:44:56 -0000 1.5 +++ TestRunTrackReconTest.java 3 Aug 2013 00:57:03 -0000 1.6 @@ -15,6 +15,7 @@
import java.util.List; import java.util.logging.Level; import java.util.logging.Logger;
+import static junit.framework.Assert.assertTrue;
import junit.framework.TestCase; import org.lcsim.event.EventHeader; import org.lcsim.event.Track;
@@ -32,7 +33,7 @@
* @author phansson */ public class TestRunTrackReconTest extends TestCase {
-
+
static final String testURLBase = "http://www.slac.stanford.edu/~phansson/files/hps_java_test/HPSTestRunv3/hps-java-1.7-SNAPSHOT-050113"; static final String testFileName = "egs_5.5gev_0.016x0_500mb_recoil_recon_1_hpsTestRunTrackingTest.slcio"; static final String testURLBaseCmp = "http://www.slac.stanford.edu/~phansson/files/hps_java_test/HPSTestRunv3/hps-java-1.7-SNAPSHOT-050113";
@@ -41,15 +42,15 @@
static final boolean saveForReference = false; static final boolean cmpHistograms = true; private final int nEvents = 5000;
-
+
public void testTrackRecon() throws Exception {
-
+
File lcioInputFile = null;
-
+
URL testURL = new URL(testURLBase + "/" + testFileName); FileCache cache = new FileCache();
- lcioInputFile = cache.getCachedFile(testURL); -
+ lcioInputFile = cache.getCachedFile(testURL); +
//Process and write out the file LCSimLoop loop = new LCSimLoop(); loop.setLCIORecordSource(lcioInputFile);
@@ -59,7 +60,7 @@
loop.add(new LCIODriver(outputFile)); loop.loop(nEvents, null); loop.dispose();
-
+
//Read LCIO back and test! LCSimLoop readLoop = new LCSimLoop(); readLoop.add(new TrackReconTestDriver());
@@ -67,32 +68,26 @@
readLoop.loop(nEvents, null); readLoop.dispose(); }
- -
+
/* * static nested class that runs the drivers for standard HPS track recon */
- - -
static class TrackReconTestDriver extends Driver { private AIDA aida = AIDA.defaultInstance(); private IAnalysisFactory af = aida.analysisFactory();
- private final IHistogram1D hntracks = aida.histogram1D("hntracks",10,0,10); - private final IHistogram1D hnstereohits = aida.histogram1D("hnstereohits",20,0,20); - private final IHistogram1D hhitpositionx = aida.histogram1D("hhitpositionx",50,0,800); - private final IHistogram1D hhitpositiony = aida.histogram1D("hhitpositiony",50,-100,100); - private final IHistogram1D hhitpositionz = aida.histogram1D("hhitpositionz",50,-100,100); - private final List<String> histograms = Arrays.asList("hntracks","hnstereohits","hhitpositionx","hhitpositiony","hhitpositionz");
+ private final IHistogram1D hntracks = aida.histogram1D("hntracks", 10, 0, 10); + private final IHistogram1D hnstereohits = aida.histogram1D("hnstereohits", 20, 0, 20); + private final IHistogram1D hhitpositionx = aida.histogram1D("hhitpositionx", 50, 0, 800); + private final IHistogram1D hhitpositiony = aida.histogram1D("hhitpositiony", 50, -100, 100); + private final IHistogram1D hhitpositionz = aida.histogram1D("hhitpositionz", 50, -100, 100); + private final List<String> histograms = Arrays.asList("hntracks", "hnstereohits", "hhitpositionx", "hhitpositiony", "hhitpositionz");
static final double alpha = 0.05; // Type-I error rate private int ntracks = 0; private int nevents = 0;
-
//Test thresholds private final double ftracks_thr = 0.1;
- -
+
@Override protected void process(EventHeader event) { super.process(event);
@@ -100,25 +95,25 @@
//System.out.printf("%s: found %d tracks\n",this.getName(),tracks.size()); hntracks.fill(tracks.size()); ntracks += tracks.size();
- for(Track track : tracks) {
+ for (Track track : tracks) {
List<TrackerHit> hitsOnTrack = track.getTrackerHits(); hnstereohits.fill(hitsOnTrack.size());
- for(TrackerHit hit : hitsOnTrack) {
+ for (TrackerHit hit : hitsOnTrack) {
double pos[] = hit.getPosition(); this.hhitpositionx.fill(pos[0]); this.hhitpositiony.fill(pos[1]); this.hhitpositionz.fill(pos[2]); } }
-
+
++nevents; } @Override protected void endOfData() { super.endOfData();
- - if(saveForReference) {
+ + if (saveForReference) {
File outputFile = new TestOutputFile(testFileName.replaceAll(".slcio", "") + "_hpsTestRunTrackingTest.aida"); try { aida.saveAs(outputFile);
@@ -126,25 +121,25 @@
Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex); } }
-
+
/* * Basic tests */
- assertTrue("Failed to find any tracks",ntracks>0); - double ftracks = ((double)ntracks)/(double)nevents; - assertTrue("Failed to reconstruct more than " + this.ftracks_thr + " of tracks/event ("+ftracks+")",ftracks>this.ftracks_thr); - assertTrue("Failed to find any stereo hits",this.hnstereohits.mean()>0.); - - -
+ assertTrue("Failed to find any tracks", ntracks > 0); + double ftracks = ((double) ntracks) / (double) nevents; + assertTrue("Failed to reconstruct more than " + this.ftracks_thr + " of tracks/event (" + ftracks + ")", ftracks > this.ftracks_thr); + assertTrue("Failed to find any stereo hits", this.hnstereohits.mean() > 0.); + + +
// IPlotter plotter = af.createPlotterFactory().create(); // plotter.createRegions(1, 3, 0); // plotter.setTitle("Nr of tracks"); // plotter.style().statisticsBoxStyle().setVisible(false); // plotter.region(0).plot(hntracks); // plotter.show();
- - if(cmpHistograms) {
+ + if (cmpHistograms) {
File aidaCmpInputFile = null;
@@ -154,12 +149,11 @@
FileCache cache = null; try { cache = new FileCache();
-aida.saveAs("testResults.aida");
} catch (IOException ex) { Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex); } try {
- aidaCmpInputFile = cache.getCachedFile(cmpURL);
+ aidaCmpInputFile = cache.getCachedFile(cmpURL);
} catch (IOException ex) { Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex); }
@@ -171,23 +165,28 @@
} catch (IOException ex) { Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex); }
- if(tree_cmp==null) {
+ if (tree_cmp == null) {
throw new RuntimeException("cannot create the ITree for the comparison file located at" + aidaCmpInputFile.getAbsolutePath()); }
- for(String histname : histograms) { - - IHistogram1D h_ref = (IHistogram1D)tree_cmp.find(histname);
+ for (String histname : histograms) { + + IHistogram1D h_ref = (IHistogram1D) tree_cmp.find(histname);
IHistogram1D h_test = aida.histogram1D(histname);
- boolean nullHypoIsRejected = CompareHistograms.instance().getTTest(alpha, h_test.mean(), h_ref.mean(), h_test.rms()*h_test.rms(), h_ref.rms()*h_ref.rms(), h_test.allEntries(), h_ref.allEntries()); - double p_value = CompareHistograms.instance().getTTestPValue(h_test.mean(), h_ref.mean(), h_test.rms()*h_test.rms(), h_ref.rms()*h_ref.rms(), h_test.allEntries(), h_ref.allEntries());
+ boolean nullHypoIsRejected = CompareHistograms.instance().getTTest(alpha, h_test.mean(), h_ref.mean(), h_test.rms() * h_test.rms(), h_ref.rms() * h_ref.rms(), h_test.allEntries(), h_ref.allEntries()); + double p_value = CompareHistograms.instance().getTTestPValue(h_test.mean(), h_ref.mean(), h_test.rms() * h_test.rms(), h_ref.rms() * h_ref.rms(), h_test.allEntries(), h_ref.allEntries()); + System.out.printf("%s: %s %s T-Test (%.1f%s C.L.) with p-value=%.3f\n", TestRunTrackReconTest.class.getName(), histname, (nullHypoIsRejected ? "FAILED" : "PASSED"), (1 - alpha) * 100, "%", p_value); + assertTrue("Failed T-Test (" + (1 - alpha) * 100 + "% C.L. p-value=" + p_value + ") comparing histogram " + histname, !nullHypoIsRejected); + double ks_p_value = CompareHistograms.getKolmogorovPValue(h_test, h_ref); + boolean ksNullHypoIsRejected = (ks_p_value < alpha); + System.out.printf("%s: %s %s Kolmogorov-Smirnov test (%.1f%s C.L.) with p-value=%.3f\n", TestRunTrackReconTest.class.getName(), histname, (ksNullHypoIsRejected ? "FAILED" : "PASSED"), (1 - alpha) * 100, "%", ks_p_value); + assertTrue("Failed Kolmogorov-Smirnov test (" + (1 - alpha) * 100 + "% C.L. p-value=" + ks_p_value + ") comparing histogram " + histname, !ksNullHypoIsRejected); +
//TODO: use a real two-sample Poisson test // boolean entriesInconsistent = CompareHistograms.instance().getTTest(alpha, h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries()); // double p_value_entries = CompareHistograms.instance().getTTestPValue(h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries(), h_test.entries(), h_ref.entries());
- System.out.printf("%s: %s %s T-Test (%.1f%s C.L.) with p-value=%.3f\n",TestRunTrackReconTest.class.getName(),histname,(nullHypoIsRejected?"FAILED":"PASSED"),(1-alpha)*100,"%",p_value);
// System.out.printf("%s: %s entries are %s (N=%d Ref=%d for T-Test w/ %.1f%s C.L.)\n",TestRunTrackReconTest.class.getName(),histname,(entriesInconsistent?"INCONSISTENT":"CONSISTENT"),h_test.entries(),h_ref.entries(),(1-alpha)*100,"%");
- assertTrue("Failed T-Test ("+ (1-alpha)*100 + "% C.L. p-value=" + p_value + ") comparing histogram " + histname,!nullHypoIsRejected);
// assertTrue("Failed T-Test ("+ (1-alpha)*100 + "% C.L. p-value=" + p_value_entries + ") on number of entries comparing histogram " + histname,!entriesInconsistent);
-
+
}
@@ -195,13 +194,7 @@
Logger.getLogger(TestRunTrackReconTest.class.getName()).log(Level.SEVERE, null, ex); } }
- - }
- - -
+ }
}
- -
}
diff -u -r1.1 -r1.2 --- CompareHistograms.java 28 Nov 2012 06:36:29 -0000 1.1 +++ CompareHistograms.java 3 Aug 2013 00:57:03 -0000 1.2 @@ -1,42 +1,42 @@
-/* - * To change this template, choose Tools | Templates - * and open the template in the editor. - */
package org.lcsim.hps.util;
+import hep.aida.IAxis; +import hep.aida.IHistogram1D;
import java.util.logging.Level; import java.util.logging.Logger; import org.apache.commons.math.MathException; import org.apache.commons.math.stat.descriptive.StatisticalSummaryValues; import org.apache.commons.math.stat.inference.TTestImpl;
+import org.apache.commons.math3.distribution.KolmogorovSmirnovDistribution;
/** *
- * Class to do various comparisons of histograms - * Singleton instance with lazy instantiation - *
+ * Class to do various comparisons of histograms Singleton instance with lazy + * instantiation + *
* @author phansson */ public class CompareHistograms {
-
+
public static CompareHistograms _instance = null; TTestImpl tTest;
+
private CompareHistograms() {
- tTest = new TTestImpl();
+ tTest = new TTestImpl();
}
-
+
public static CompareHistograms instance() {
- if(_instance==null) {
+ if (_instance == null) {
_instance = new CompareHistograms(); } return _instance; }
-
+
public double getTTestPValue(double m1, double m2, double v1, double v2, int n1, int n2) {
- StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1,n1, 1.,0., 0.); - StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2,n2, 1.,0., 0.);
+ StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1, n1, 1., 0., 0.); + StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2, n2, 1., 0., 0.);
- double p_value = -1;
+ double p_value = -1;
try { p_value = tTest.tTest(stat1, stat2); } catch (IllegalArgumentException ex) {
@@ -45,12 +45,12 @@
Logger.getLogger(CompareHistograms.class.getName()).log(Level.SEVERE, null, ex); } return p_value;
-
+
}
- public boolean getTTest(double alpha,double m1, double m2, double v1, double v2, int n1, int n2) { - StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1,n1, 1.,0., 0.); - StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2,n2, 1.,0., 0.);
+ public boolean getTTest(double alpha, double m1, double m2, double v1, double v2, int n1, int n2) { + StatisticalSummaryValues stat1 = new StatisticalSummaryValues(m1, v1, n1, 1., 0., 0.); + StatisticalSummaryValues stat2 = new StatisticalSummaryValues(m2, v2, n2, 1., 0., 0.);
boolean nullHypoIsRejected = false; try {
@@ -62,8 +62,26 @@
} return nullHypoIsRejected; }
- -
-
+ public static double getKolmogorovPValue(IHistogram1D reference, IHistogram1D test) { + double integralMax = 0.0; + double refIntegral = reference.binHeight(IAxis.UNDERFLOW_BIN); + double testIntegral = test.binHeight(IAxis.UNDERFLOW_BIN); + double integralDiff = Math.abs(refIntegral / reference.allEntries() - testIntegral / test.allEntries()); + if (integralDiff > integralMax) { + integralMax = integralDiff; + } + for (int i = 0; i < reference.axis().bins(); i++) { + refIntegral += reference.binHeight(i); + testIntegral += test.binHeight(i); + + integralDiff = Math.abs(refIntegral / reference.allEntries() - testIntegral / test.allEntries()); + if (integralDiff > integralMax) { + integralMax = integralDiff; + } + } + int n = (int) Math.ceil(Math.sqrt((reference.allEntries() * test.allEntries()) / (reference.allEntries() + test.allEntries()))); + KolmogorovSmirnovDistribution dist = new KolmogorovSmirnovDistribution(n); + return 1.0 - dist.cdf(integralMax); + }
}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1