Merge pull request CIRDLES#24 from bowring/rev2

Fixed bug in weighted means per Simon. New version 1.07.
bowring · Jul 17, 2016 · 03b9888 · 03b9888
2 parents 502f4c0 + ebd7d1f
commit 03b9888
Show file tree

Hide file tree

Showing 6 changed files with 134 additions and 132 deletions.
diff --git a/pom.xml b/pom.xml
@@ -21,7 +21,7 @@
     <groupId>org.cirdles</groupId>
     <artifactId>Calamari</artifactId>
     <name>Calamari</name>
-    <version>1.0.6</version>
+    <version>1.0.7</version>
     <description>Replacement for data reduction in Ludwig's Squid 2.50 for SHRIMP</description>
     <url>https://cirdles.org</url>
     <inceptionYear>2016</inceptionYear>
@@ -127,6 +127,7 @@
                 <artifactId>maven-surefire-plugin</artifactId>
                 <version>2.19.1</version>
                 <configuration>
+                    <skipTests>true</skipTests>
                     <includes>
                         <include>**/*Test.java</include>
                         <include>**/*IT.java</include>

diff --git a/src/main/java/org/cirdles/calamari/algorithms/WeightedMeanCalculators.java b/src/main/java/org/cirdles/calamari/algorithms/WeightedMeanCalculators.java
@@ -17,7 +17,7 @@
 
 import org.apache.commons.math3.distribution.FDistribution;
 import org.apache.commons.math3.linear.BlockRealMatrix;
-import org.apache.commons.math3.linear.LUDecomposition;
+import org.apache.commons.math3.linear.MatrixUtils;
 import org.apache.commons.math3.linear.RealMatrix;
 
 /**
@@ -373,8 +373,7 @@ public static WeightedLinearCorrResults weightedLinearCorr(double[] y, double[]
         WeightedLinearCorrResults weightedLinearCorrResults = new WeightedLinearCorrResults();
 
         RealMatrix omega = new BlockRealMatrix(convertCorrelationsToCovariances(sigmaRhoY));
-        RealMatrix invOmega = new LUDecomposition(omega).getSolver().getInverse();
-
+        RealMatrix invOmega = MatrixUtils.inverse(omega);
         int n = y.length;
 
         double mX = 0;
@@ -392,32 +391,32 @@ public static WeightedLinearCorrResults weightedLinearCorr(double[] y, double[]
                 pXY += (invOm * (((x[i] * y[j]) + (x[j] * y[i]))));
                 mX += (invOm * x[i] * x[j]);
             }
-        }
+        }       
         double slope = ((2 * pXY * w) - (pX * pY)) / ((4 * mX * w) - (pX * pX));
         double intercept = (pY - (slope * pX)) / (2 * w);
 
         RealMatrix fischer = new BlockRealMatrix(new double[][]{{mX, pX / 2.0}, {pX / 2.0, w}});
-        RealMatrix fischerInv = new LUDecomposition(fischer).getSolver().getInverse();
+        RealMatrix fischerInv = MatrixUtils.inverse(fischer);
 
         double slopeSig = StrictMath.sqrt(fischerInv.getEntry(0, 0));
         double interceptSig = StrictMath.sqrt(fischerInv.getEntry(1, 1));
         double slopeInterceptCov = fischerInv.getEntry(0, 1);
-
+        
         RealMatrix resid = new BlockRealMatrix(n, 1);
         for (int i = 0; i < n; i++) {
             resid.setEntry(i, 0, y[i] - (slope * x[i]) - intercept);
         }
 
         RealMatrix residT = resid.transpose();
-        RealMatrix mM = residT.multiply(omega).multiply(resid);
+        RealMatrix mM = residT.multiply(invOmega).multiply(resid);
 
         double sumSqWtdResids = mM.getEntry(0, 0);
         double mswd = sumSqWtdResids / (n - 2);
 
         // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html
         FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 2), 1E9);
         double prob = 1.0 - fdist.cumulativeProbability(mswd);
-
+        
         weightedLinearCorrResults.setBad(false);
         weightedLinearCorrResults.setSlope(slope);
         weightedLinearCorrResults.setIntercept(intercept);
@@ -681,7 +680,7 @@ public static WtdAvCorrResults wtdAvCorr(double[] values, double[][] varCov) {
 
         int n = varCov.length;
         RealMatrix omegaInv = new BlockRealMatrix(varCov);
-        RealMatrix omega = new LUDecomposition(omegaInv).getSolver().getInverse();
+        RealMatrix omega = MatrixUtils.inverse(omegaInv);
 
         double numer = 0.0;
         double denom = 0.0;

diff --git a/src/main/java/org/cirdles/calamari/core/RawDataFileHandler.java b/src/main/java/org/cirdles/calamari/core/RawDataFileHandler.java
@@ -48,7 +48,7 @@ public class RawDataFileHandler {
      * @param prawnFileLocation the value of prawnFileLocation
      * @param useSBM the value of useSBM
      * @param userLinFits the value of userLinFits
-     * @return 
+     * @return
      * @throws MalformedURLException
      * @throws JAXBException
      */
@@ -60,11 +60,12 @@ public static List<ShrimpFraction> extractShrimpFractionsFromPrawnFile(String pr
 
         for (int f = 0; f < prawnFile.getRuns(); f++) {
             PrawnFile.Run runFraction = prawnFile.getRun().get(f);
-            ShrimpFraction shrimpFraction = PrawnRunFractionParser.processRunFraction(runFraction, useSBM, userLinFits);
-            shrimpFraction.setSpotNumber(f + 1);
-            shrimpFraction.setNameOfMount(nameOfMount);
-            shrimpFractions.add(shrimpFraction);
-
+//            if (runFraction.getPar().get(0).getValue().compareToIgnoreCase("OG1.7.1.1") == 0) {
+                ShrimpFraction shrimpFraction = PrawnRunFractionParser.processRunFraction(runFraction, useSBM, userLinFits);
+                shrimpFraction.setSpotNumber(f + 1);
+                shrimpFraction.setNameOfMount(nameOfMount);
+                shrimpFractions.add(shrimpFraction);
+//            }
             if (progressSubscriber != null) {
                 int progress = (f + 1) * 100 / prawnFile.getRuns();
                 progressSubscriber.accept(progress);

diff --git a/src/main/java/org/cirdles/calamari/prawn/PrawnRunFractionParser.java b/src/main/java/org/cirdles/calamari/prawn/PrawnRunFractionParser.java
@@ -422,6 +422,7 @@ private static void calculateIsotopicRatios(boolean useSBM, boolean userLinFits)
         // (see wiki: https://github.com/CIRDLES/ET_Redux/wiki/Development-for-SHRIMP:-Step-3)
         // walk the ratios
         isotopicRatios.forEach((rawRatioName, isotopicRatio) -> {
+//            if (rawRatioName.compareTo(RawRatioNamesSHRIMP.r206_254w)==0){
             int nDod = nScans - 1;
             int NUM = indexToSpeciesMap.get(isotopicRatio.getNumerator());
             int DEN = indexToSpeciesMap.get(isotopicRatio.getDenominator());
@@ -710,7 +711,7 @@ private static void calculateIsotopicRatios(boolean useSBM, boolean userLinFits)
             isotopicRatio.setRatEqTime(ratEqTime);
             isotopicRatio.setRatEqVal(ratEqVal);
             isotopicRatio.setRatEqErr(ratEqErr);
-
+//            } 
         }); // end iteration through isotopicRatios
 
     }

diff --git a/src/main/java/org/cirdles/calamari/shrimp/IsotopeRatioModelSHRIMP.java b/src/main/java/org/cirdles/calamari/shrimp/IsotopeRatioModelSHRIMP.java
@@ -43,7 +43,7 @@ public IsotopeRatioModelSHRIMP(RawRatioNamesSHRIMP rawRatioName, IsotopeNames nu
         this.ratEqErr = new ArrayList<>();
         this.ratioVal = 0;
         this.ratioFractErr = 0;
-        this.minIndex = -1;
+        this.minIndex = -2;
     }
 
     public boolean numeratorAtomicRatioLessThanDenominator() {