added some error messages for data input format problems in model/data.d

stschiff · Jan 23, 2015 · 26733ca · 26733ca
1 parent ca88bd3
commit 26733ca
Show file tree

Hide file tree

Showing 4 changed files with 23 additions and 17 deletions.
diff --git a/Makefile b/Makefile
@@ -28,7 +28,7 @@ testcoverage : model/*.d unittest.d powell.d brent.d maximization_step.d expecta
 	build/unittest
 	mv *.lst code_coverage/
 
-unittest : model/*.d unittest.d powell.d brent.d maximization_step.d expectation_step.d amoeba.d logger.d branchlength.d
+unittest : model/*.d unittest.d powell.d brent.d maximization_step.d expectation_step.d logger.d branchlength.d
 	dmd -unittest ${GSL} -odbuild -ofbuild/unittest $^
 	build/unittest
 

diff --git a/maximization_step.d b/maximization_step.d
@@ -241,7 +241,7 @@ unittest {
   auto minFunc = new MinFunc(expectationResultVec, expectationResultMat, params, timeSegmentPattern, false, false);
   auto rho = 0.001;
   auto x = minFunc.getXfromLambdaVec(lambdaVec);
-  x ~= minFunc.getXfromRecombinationRate(rho);
+  x ~= log(rho);
   auto lambdaFromX = minFunc.getLambdaVecFromX(x);
   auto rhoFromX = minFunc.getRecombinationRateFromX(x);
   foreach(i; 0 .. lambdaVec.length)
@@ -250,7 +250,7 @@ unittest {
 
   minFunc = new MinFunc(expectationResultVec, expectationResultMat, params, timeSegmentPattern, true, false);
   x = minFunc.getXfromLambdaVec(lambdaVec);
-  x ~= minFunc.getXfromRecombinationRate(rho);
+  x ~= log(rho);
   lambdaFromX = minFunc.getLambdaVecFromXfixedPop(x);
   rhoFromX = minFunc.getRecombinationRateFromX(x);
   foreach(i; 0 .. lambdaVec.length)

diff --git a/model/data.d b/model/data.d
@@ -114,26 +114,27 @@ unittest {
   assert(normalizeAlleleString("TC") == "01");
 }
 
-void checkDataLine(const char[] line) {
-  auto r = regex(r"^\w+\s\d+\s\d+(\s[ACTG01\?,]+){0,1}$");
-  enforce(match(line, r));
-}
-
-unittest {
-  assertThrown(checkDataLine("1 20 5 AACC,AACA 2.44"));
-  assertNotThrown(checkDataLine("1 20 5 AACC"));
-  assertNotThrown(checkDataLine("4 5 2"));
-  assertNotThrown(checkDataLine("1 10 5 ACC"));
-  assertThrown(checkDataLine("1 20 5 AGGSSXX"));
-}
+// void checkDataLine(const char[] line) {
+//   auto r = regex(r"^\w+\s\d+\s\d+(\s[ACTG01\?,]+){0,1}$");
+//   enforce(match(line, r));
+// }
+//
+// unittest {
+//   assertThrown(checkDataLine("1 20 5 AACC,AACA 2.44"));
+//   assertNotThrown(checkDataLine("1 20 5 AACC"));
+//   assertNotThrown(checkDataLine("4 5 2"));
+//   assertNotThrown(checkDataLine("1 10 5 ACC"));
+//   assertThrown(checkDataLine("1 20 5 AGGSSXX"));
+// }
 
 size_t getNrHaplotypesFromFile(string filename) {
   auto file = File(filename, "r");
   scope(exit) file.close();
   auto line = file.readln();
   line = line.strip();
-  checkDataLine(line);
+  // checkDataLine(line);
   auto fields = line.strip().split();
+  enforce(fields.length > 2, "illegal input file");
   if(fields.length < 4)
     return 2;
   else {
@@ -178,11 +179,17 @@ SegSite_t[] readSegSites(string filename, bool directedEmissions, size_t[] indic
 
   auto f = File(filename, "r");
   long lastPos = -1;
+  auto chrom = "";
   foreach(line; f.byLine()) {
     // checkDataLine(line.strip());
     auto fields = line.strip().split();
+    if(chrom == "")
+        chrom = fields[0].idup;
+    else
+        enforce(chrom == fields[0], "chromosomes must all be the same within one file (sorry)");
     auto pos = to!size_t(fields[1]);
     auto nrCalledSites = to!size_t(fields[2]);
+    enforce(nrCalledSites > 0, "nr of called sites (3rd column in input file) must be always > 0");
     if(lastPos == -1) {
       lastPos = pos - nrCalledSites;
     }

diff --git a/unittest.d b/unittest.d
@@ -38,7 +38,6 @@ import powell;
 import brent;
 import maximization_step;
 import expectation_step;
-import amoeba;
 import branchlength;
 
 void main() {