diff --git a/README.md b/README.md
index 2d9ddb47f6d50e8e2b7a255dca3284813f63868d..3b59162aa963268ce274afab8ec29166e517fe30 100644
--- a/README.md
+++ b/README.md
@@ -1,41 +1,43 @@
# REQ
-_REQ_ is a command line program written in [Java](https://docs.oracle.com/javase/8/docs/technotes/guides/language/index.html) that allows estimating the rate of elementary quartets (REQ) for each internal branch of a phylogenetic tree from a distance matrix, as described by [Guénoche and Garreta (2001)](https://doi.org/10.1007/3-540-45727-5_5).
+[_REQ_](https://research.pasteur.fr/fr/tool/r%CE%B5q-assessing-branch-supports-o%C6%92-a-distance-based-phylogenetic-tree-with-the-rate-o%C6%92-elementary-quartets/) is a command line program written in [Java](https://docs.oracle.com/javase/8/docs/technotes/guides/language/index.html) that estimates the rate of elementary quartets (REQ) for each internal branch of an unrooted binary phylogenetic tree from a distance matrix, as described by [Guénoche and Garreta (2001)](https://doi.org/10.1007/3-540-45727-5_5).
## Compilation and execution
-The source code of _REQ_ is inside the _src_ directory and could be compiled and executed in two different ways.
+The source code of _REQ_ is inside the _src_ directory and can be compiled and executed in two different ways.
#### Building an executable jar file
-On computers with [Oracle JDK](http://www.oracle.com/technetwork/java/javase/downloads/index.html) (6 or higher) installed, a Java executable jar file could be created. In a command-line window, go to the _src_ directory and type:
+On computers with [Oracle JDK](http://www.oracle.com/technetwork/java/javase/downloads/index.html) (6 or higher) installed, a Java executable jar file can be created.
+In a command-line window, go to the _src_ directory and type:
```bash
javac REQ.java
echo Main-Class: REQ > MANIFEST.MF
jar -cmvf MANIFEST.MF REQ.jar REQ.class
rm MANIFEST.MF REQ.class
```
-This will create the executable jar file `REQ.jar` that could be launched with the following command line model:
+This will create the executable jar file `REQ.jar` that can be run with the following command line model:
```bash
java -jar REQ.jar [files]
```
#### Building a native executable
-On computers with [GraalVM](hhttps://www.graalvm.org/downloads/) installed, a native executable could also be built. In a command-line window, go to the _src_ directory, and type:
+On computers with [GraalVM](hhttps://www.graalvm.org/downloads/) installed, a native executable can also be built.
+In a command-line window, go to the _src_ directory, and type:
```bash
javac REQ.java
native-image REQ REQ
rm REQ.class
```
-This will create the native executable `REQ` that could be launched with the following command line model:
+This will create the native executable `REQ` that can be launched with the following command line model:
```bash
./REQ [files]
```
## Usage
-Launch _REQ_ without option to read the following documentation:
+Run _REQ_ without option to read the following documentation:
```
USAGE: REQ <dfile> <tfile> <outfile> [-v]
@@ -44,7 +46,7 @@ Launch _REQ_ without option to read the following documentation:
square format
<tfile> unrooted binary phylogenetic tree file with no confidence
value at branches in NEWICK format
- <outfile> outfile name
+ <outfile> output file name
-v verbose mode
```
@@ -58,7 +60,7 @@ The following command line writes into the file _tree.req.t_ the phylogenetic tr
```bash
REQ example/matrix.d example/tree.t tree.req.t -v
```
-Because the option `-v` is set, the verbose mode will output the tree topology in NEWICK format, the list of leaf names, and, for each internal branch, the leaf quadripartition following by the rate of elementary quartets _Re_:
+Because the option `-v` is set, the verbose mode will output the tree topology in NEWICK format, the list of leaf names, and, for each internal branch, the leaf quadripartition followed by the rate of elementary quartets _Re_:
```
# (((((((17,18),16),((20,21),19)),(((4,((6,7),5)),(((2,3),1),0)),(8,9))),(10,11)),(12,13)),14,15);
0: P07_621_SLS
@@ -83,25 +85,25 @@ Because the option `-v` is set, the verbose mode will output the tree topology i
19: P03_594_BU1_SLS
20: P05_599_BU1_SLS
21: P03_592_BU1_SLS
-[17][18][16][21,20,19,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0] Re=1.000 (19/19)
-[18,17][16][19,21,20][15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0] Re=0.677 (65/96)
-[20][21][19][18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0] Re=0.474 (9/19)
-[21,20][19][16,18,17][15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0] Re=1.000 (96/96)
-[16,18,17][19,21,20][9,8,0,1,3,2,5,7,6,4][15,14,13,12,11,10] Re=0.856 (462/540)
-[6][7][5][21,20,19,18,17,16,15,14,13,12,11,10,9,8,4,3,2,1,0] Re=1.000 (19/19)
-[7,6][5][4][21,20,19,18,17,16,15,14,13,12,11,10,9,8,3,2,1,0] Re=0.944 (34/36)
-[4][5,7,6][0,1,3,2][21,20,19,18,17,16,15,14,13,12,11,10,9,8] Re=1.000 (168/168)
-[2][3][1][21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,0] Re=1.000 (19/19)
-[3,2][1][0][21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4] Re=1.000 (36/36)
-[1,3,2][0][5,7,6,4][21,20,19,18,17,16,15,14,13,12,11,10,9,8] Re=1.000 (168/168)
-[5,7,6,4][0,1,3,2][9,8][21,20,19,18,17,16,15,14,13,12,11,10] Re=1.000 (384/384)
-[8][9][0,1,3,2,5,7,6,4][21,20,19,18,17,16,15,14,13,12,11,10] Re=1.000 (96/96)
-[0,1,3,2,5,7,6,4][9,8][19,21,20,16,18,17][15,14,13,12,11,10] Re=0.453 (261/576)
-[19,21,20,16,18,17][9,8,0,1,3,2,5,7,6,4][11,10][15,14,13,12] Re=0.487 (234/480)
-[10][11][9,8,0,1,3,2,5,7,6,4,19,21,20,16,18,17][15,14,13,12] Re=1.000 (64/64)
-[9,8,0,1,3,2,5,7,6,4,19,21,20,16,18,17][11,10][13,12][15,14] Re=0.594 (76/128)
-[12][13][11,10,9,8,0,1,3,2,5,7,6,4,19,21,20,16,18,17][15,14] Re=1.000 (36/36)
-[11,10,9,8,0,1,3,2,5,7,6,4,19,21,20,16,18,17][13,12][14][15] Re=0.500 (18/36)
+[17][18][16][0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,19,20,21] Re=1.000 (19/19)
+[17,18][16][19,20,21][0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] Re=0.677 (65/96)
+[20][21][19][0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18] Re=0.474 (9/19)
+[20,21][19][16,17,18][0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] Re=1.000 (96/96)
+[16,17,18][19,20,21][0,1,2,3,4,5,6,7,8,9][10,11,12,13,14,15] Re=0.856 (462/540)
+[6][7][5][0,1,2,3,4,8,9,10,11,12,13,14,15,16,17,18,19,20,21] Re=1.000 (19/19)
+[6,7][5][4][0,1,2,3,8,9,10,11,12,13,14,15,16,17,18,19,20,21] Re=0.944 (34/36)
+[4][5,6,7][0,1,2,3][8,9,10,11,12,13,14,15,16,17,18,19,20,21] Re=1.000 (168/168)
+[2][3][1][0,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21] Re=1.000 (19/19)
+[2,3][1][0][4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21] Re=1.000 (36/36)
+[1,2,3][0][4,5,6,7][8,9,10,11,12,13,14,15,16,17,18,19,20,21] Re=1.000 (168/168)
+[4,5,6,7][0,1,2,3][8,9][10,11,12,13,14,15,16,17,18,19,20,21] Re=1.000 (384/384)
+[8][9][0,1,2,3,4,5,6,7][10,11,12,13,14,15,16,17,18,19,20,21] Re=1.000 (96/96)
+[0,1,2,3,4,5,6,7][8,9][16,17,18,19,20,21][10,11,12,13,14,15] Re=0.453 (261/576)
+[16,17,18,19,20,21][0,1,2,3,4,5,6,7,8,9][10,11][12,13,14,15] Re=0.488 (234/480)
+[10][11][0,1,2,3,4,5,6,7,8,9,16,17,18,19,20,21][12,13,14,15] Re=1.000 (64/64)
+[0,1,2,3,4,5,6,7,8,9,16,17,18,19,20,21][10,11][12,13][14,15] Re=0.594 (76/128)
+[12][13][0,1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,20,21][14,15] Re=1.000 (36/36)
+[0,1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,20,21][12,13][14][15] Re=0.500 (18/36)
output tree written into tree.req.t
```
diff --git a/example/tree.req.t b/example/tree.req.t
index 559eebf15694cd1e7434e1cab74b1ebbd21d8b76..a28116fdc0b807e0873cff55a3274e140a063c3a 100644
--- a/example/tree.req.t
+++ b/example/tree.req.t
@@ -1 +1 @@
-(((((((P10_703_BU2_PER:0.000157849066,P09_704_BU2_PER:0.000209200934)1.000:0.000483122851,P05_622_BU1_SLS:0.001116966803)0.677:0.000018368078,((P05_599_BU1_SLS:0.000159586258,P03_592_BU1_SLS:0.000169023742)0.474:0.000015792325,P03_594_BU1_SLS:0.000128336869)1.000:0.000462006451)0.856:0.000081437521,(((P11_702_BU2_PER:0.000615002481,((P04_601_BU1_SLS:0.000117392803,P04_559_BU1_SLS:0.000164857197)1.000:0.000925936953,E01_615_SLS:0.000884566733)0.944:0.000161626490)1.000:0.008730158578,(((P08_701_BU2_PER:0.001390856250,P01_617_BU1_SLS:0.001385183750)1.000:0.001280565582,1006PHL:0.002499367591)1.000:0.011475862713,P07_621_SLS:0.014634622105)1.000:0.015451073655)1.000:0.008164485602,(P05_600_BU1_SLS:0.000172301484,P05_598_BU1_SLS:0.000210108516)1.000:0.000603348539)0.453:0.000010417106)0.488:0.000045431349,(P04_603_BU1_SLS:0.000255598014,P04_602_BU1_SLS:0.000179821986)1.000:0.000552675249)0.594:0.000029208891,(P06_032_BU1_SLS:0.000173131711,P06_023_BU1_SLS:0.000108148289)1.000:0.000554756800)0.500:0.000022656931,P02_783_BU1_SLS:0.000592772460,P12_579_STR:0.000673747540);
\ No newline at end of file
+(((((((P10_703_BU2_PER:0.000157849066,P09_704_BU2_PER:0.000209200934)1.000:0.000483122851,P05_622_BU1_SLS:0.001116966803)0.677:0.000018368078,((P05_599_BU1_SLS:0.000159586258,P03_592_BU1_SLS:0.000169023742)0.474:0.000015792325,P03_594_BU1_SLS:0.000128336869)1.000:0.000462006451)0.856:0.000081437521,(((P11_702_BU2_PER:0.000615002481,((P04_601_BU1_SLS:0.000117392803,P04_559_BU1_SLS:0.000164857197)1.000:0.000925936953,E01_615_SLS:0.000884566733)0.944:0.000161626490)1.000:0.008730158578,(((P08_701_BU2_PER:0.001390856250,P01_617_BU1_SLS:0.001385183750)1.000:0.001280565582,1006PHL:0.002499367591)1.000:0.011475862713,P07_621_SLS:0.014634622105)1.000:0.015451073655)1.000:0.008164485602,(P05_600_BU1_SLS:0.000172301484,P05_598_BU1_SLS:0.000210108516)1.000:0.000603348539)0.453:0.000010417106)0.488:0.000045431349,(P04_603_BU1_SLS:0.000255598014,P04_602_BU1_SLS:0.000179821986)1.000:0.000552675249)0.594:0.000029208891,(P06_032_BU1_SLS:0.000173131711,P06_023_BU1_SLS:0.000108148289)1.000:0.000554756800)0.500:0.000022656931,P02_783_BU1_SLS:0.000592772460,P12_579_STR:0.000673747540);
diff --git a/src/REQ.java b/src/REQ.java
index 672cca65a70e7a6083ea15f56bad42b2eb3ae78e..4f85213796ddbe4d064dfea17fbd4e236391a2f2 100644
--- a/src/REQ.java
+++ b/src/REQ.java
@@ -1,31 +1,37 @@
/*
- ####################################################################
- REQ: estimating the rate of elementary quartets (REQ) for each
- branch of a phylogenetic tree from a distance matrix
+ ########################################################################################################
- Copyright (C) 2017,2018,2019 Alexis Criscuolo
-
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
+ REQ: assessing branch supports of a distance-based phylogenetic tree with rates of elementary quartets
+
+ Copyright (C) 2017,2018,2019,2020 Institut Pasteur
+
+ This program is free software: you can redistribute it and/or modify it under the terms of the GNU
+ General Public License as published by the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
- This program is distributed in the hope that it will be useful, but
- WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- General Public License for more details.
+ This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
- Contact:
- Institut Pasteur
- Bioinformatics and Biostatistics Hub
- C3BI, USR 3756 IP CNRS
- Paris, FRANCE
-
- alexis.criscuolo@pasteur.fr
- ####################################################################
+ You should have received a copy of the GNU General Public License along with this program. If not, see
+ <http://www.gnu.org/licenses/>.
+
+ Contact:
+ Alexis Criscuolo alexis.criscuolo@pasteur.fr
+ Genome Informatics & Phylogenetics (GIPhy) giphy.pasteur.fr
+ Bioinformatics and Biostatistics Hub research.pasteur.fr/team/hub-giphy
+ USR 3756 IP CNRS research.pasteur.fr/team/bioinformatics-and-biostatistics-hub
+ Dpt. Biologie Computationnelle research.pasteur.fr/department/computational-biology
+ Institut Pasteur, Paris, FRANCE research.pasteur.fr
+
+ ########################################################################################################
+*/
+/*
+ ########################################################################################################
+ ### v2.0.200528ac
+ + verify that the input tree is unrooted and fully binary
+ + faster running times by minimizing String conversion
+ ########################################################################################################
*/
import java.io.*;
@@ -35,7 +41,8 @@ import java.text.*;
public class REQ {
// constants
- final static String VERSION = "v1.3.190304ac";
+ final static String VERSION = "v2.0.200528ac";
+ final static int MAX = 32760;
final static int INF = Integer.MAX_VALUE;
// io
@@ -43,23 +50,25 @@ public class REQ {
static BufferedWriter out;
static File dfile, tfile;
static boolean verbose;
- static NumberFormat df;
// data
static int n;
static ArrayList<String> lbl;
static float[][] dm;
static StringBuilder nwk, tr, tro;
+ static short[] tip, totip;
+ static ArrayList<Double> are;
// stuffs
- static int i, j, m, u, v, w, x, y, last, sup;
- static double re, ab, ac, bc, abcd, up, dn;
- static short si;
+ static short s, as, ae, bs, be, cs, ce;
+ static int d, td, i, j, m, u, v, w, x, y, last, sup;
+ static long up;
+ static double re, ab, ad, bd, abcd, dn;
+ static char ch;
static String line;
static String[] split;
- static short[] sta, stb, stc, std, sa;
- static TreeSet<Short> ts, stdts;
- static ArrayList<Double> are;
+ static short[] sta, stb, stc;
+ static BitSet std;
public static void main(String[] args) throws IOException {
@@ -74,119 +83,147 @@ public class REQ {
System.out.println(" square format");
System.out.println(" <tfile> unrooted binary phylogenetic tree file with no confidence");
System.out.println(" value at branches in NEWICK format");
- System.out.println(" <outfile> outfile name");
+ System.out.println(" <outfile> output file name");
System.out.println(" -v verbose mode");
System.out.println("");
System.exit(0);
}
//### init ###########################################################################################################################################################################
- df = NumberFormat.getNumberInstance(Locale.US); df.setGroupingUsed(false); df.setMaximumFractionDigits(3); df.setMinimumFractionDigits(3);
verbose = ( (args.length > 3) && args[3].equals("-v") ) ? true : false;
if ( verbose ) System.out.println("REQ " + VERSION);
- if ( ! (dfile = new File(args[0])).exists() ) { System.err.println("distance matrix file does not exist: " + args[0]); System.exit(1); }
- if ( ! (tfile = new File(args[1])).exists() ) { System.err.println("tree file does not exist: " + args[1]); System.exit(1); }
+ if ( ! (dfile = new File(args[0])).exists() ) { System.err.println("distance matrix file does not exist: " + args[0]); System.exit(1); }
+ if ( ! (tfile = new File(args[1])).exists() ) { System.err.println("tree file does not exist: " + args[1]); System.exit(1); }
//### reading distance matrix dm #####################################################################################################################################################
in = new BufferedReader(new FileReader(dfile));
- while ( true ) try { if ( (line=in.readLine().trim()).length() != 0 ) break; } catch ( NullPointerException e ) { System.out.println("matrix file is empty"); System.exit(1); }
- try { n = Integer.parseInt(line); } catch ( NumberFormatException e ) { System.out.println("distance matrix file is incorrectly formatted: " + args[0]); System.exit(1); }
- if ( n > 32760 ) { System.out.println("too many taxa (>32760)"); System.exit(1); }
+ while ( true ) try { if ( (line=in.readLine().trim()).length() != 0 ) break; } catch ( NullPointerException e ) { System.err.println("matrix file is empty"); System.exit(1); }
+ try { if ( (n=Integer.parseInt(line)) > MAX ) { System.err.println("too many taxa (>" + MAX + ")"); System.exit(1); } }
+ catch ( NumberFormatException e ) { System.err.println("distance matrix file is incorrectly formatted: " + line); System.exit(1); }
lbl = new ArrayList<String>(n); dm = new float[n][]; i = -1;
- while ( true ) {
- try { line = in.readLine().trim(); } catch ( NullPointerException e ) { in.close(); break; }
- split = line.split("\\s+"); lbl.add(split[0]); dm[++i] = new float[i]; j = 0; while ( ++j <= i ) dm[i][--j] = Float.parseFloat(split[++j]);
- }
- if ( ++i != n ) { System.out.println("distance matrix file is incorrectly formatted: " + args[0]); System.exit(1); }
+ try {
+ while ( ++i < n ) {
+ if ( (line=in.readLine().trim()).length() == 0 ) { --i; continue; }
+ split = line.split("\\s+"); lbl.add(split[j=0]); dm[i] = new float[i]; while ( ++j <= i ) dm[i][--j] = (float) Double.parseDouble(split[++j]);
+ }
+ in.close();
+ } catch ( NullPointerException e ) { System.err.println("distance matrix file is incorrectly formatted"); System.exit(1); }
//### reading phylogenetic tree nwk ##################################################################################################################################################
nwk = new StringBuilder(""); in = new BufferedReader(new FileReader(tfile));
- while ( true ) try { nwk = nwk.append(in.readLine().trim()); } catch ( NullPointerException e ) { in.close(); break; }
+ try { while ( true ) nwk = nwk.append(in.readLine().trim()); } catch ( NullPointerException e ) { in.close(); }
tr = new StringBuilder(nwk.toString());
//### discarding branch lengths from tr ##############################################################################################################################################
- //u = -1; while ( (u=tr.indexOf(":", ++u)) != -1 ) { x = ( (x=tr.indexOf(",", u)) != -1 ) ? x : INF; y = ( (y=tr.indexOf(")", u)) != -1 ) ? y : INF; tr = tr.delete(u, (x<y)?x:y); }
u = -1 ; while ( (u=tr.indexOf(":", ++u)) != -1 ) tr = tr.delete(u, (x=((x=tr.indexOf(",",u))!=-1)?x:INF) < (y=((y=tr.indexOf(")",u))!=-1)?y:INF) ? x : y);
//### encoding lbl within tr #########################################################################################################################################################
- u = -1;
- while ( (u = (x=((x=tr.indexOf(",",u))!=-1)?x:INF) < (y=((y=tr.indexOf("(",u))!=-1)?y:INF) ? x : y) < INF ) {
- if ( tr.charAt(++u) == '(' ) continue;
- tr = tr.replace(u, (v = (x=((x=tr.indexOf(",",u))!=-1)?x:INF) < (y=((y=tr.indexOf(")",u))!=-1)?y:INF) ? x : y), String.valueOf(j=lbl.indexOf(tr.substring(u, v))));
- if ( j == -1 ) { System.out.println(tr.substring(u, v) + " is not inside distance matrix"); System.exit(1); }
- }
- if ( verbose ) System.out.println("# " + tr.toString());
-
+ tip = new short[n]; i = n; //## NOTE: tip[i] = ith tax id in tr
+ d = 0; //## NOTE: d = no. closing parentheses in tr
+ u = tr.length(); v = --u;
+ while ( --u >= 0 )
+ switch ( tr.charAt(u) ) {
+ case '(':
+ if ( tr.charAt(++u) != '(' ) { tr = tr.replace(u, v, String.valueOf(j=lbl.indexOf(tr.substring(u, v))));
+ if ( j != -1 ) tip[--i] = (short) j; else { System.err.println(tr.substring(u, v) + " is not inside distance matrix"); System.exit(1); } }
+ --u; break;
+ case ',':
+ if ( tr.charAt(++u) != '(' ) { tr = tr.replace(u, v, String.valueOf(j=lbl.indexOf(tr.substring(u, v))));
+ if ( j != -1 ) tip[--i] = (short) j; else { System.err.println(tr.substring(u, v) + " is not inside distance matrix"); System.exit(1); } }
+ v = --u; break;
+ case ')': ++d; v = u; break;
+ }
+ if ( n != d + 2 ) { System.err.println(( n == d + 1 ) ? "input tree is rooted" : "input tree is not binary"); System.exit(1); }
+ if ( verbose ) { System.out.println("# " + tr.toString()); i = -1; while ( ++i < n ) System.out.println(i + ": " + lbl.get(i)); }
+
//### estimating the rate of EQ re for each internal branch e ########################################################################################################################
- tro = new StringBuilder(tr.toString());
+ tro = new StringBuilder(tr.toString());
//# first tr rooting: ((ST1),(ST2),(ST3)); => (((ST1),(ST2)),(ST3)); ###############################################################
//# | |
//# sup last
tr = tr.insert(0, '('); sup = apc(tr, tr.lastIndexOf(")")); tr = tr.insert(sup, ')'); //# NOTE: closing parenthesis at index 'sup' should not be considered for REQ calculations
- last = tr.lastIndexOf(")");
- if ( verbose ) { i = -1; while ( ++i < n ) System.out.println(i + ": " + lbl.get(i)); /*System.out.println("# " + tr.toString());*/ }
-
- ts = new TreeSet<Short>(); si = 0; --si; while ( ++si < n ) ts.add(new Short(si)); are = new ArrayList<Double>(); u = -1;
+ last = tr.lastIndexOf(")");
+ //# NOTE: totip[x] = id in tip[] of the taxon at or right after position x in tr //# tr = "( ( 1 , 3 ) , ( 0 , 2 ) ) ;" tip = [1,3,0,2]
+ totip = new short[x=tr.length()]; //# totip = [0,0,0,1,1,2,2,2,2,3,4,4,4,4]
+ s = (short) n; while ( --x >= 0 ) { switch ( tr.charAt(x) ) { case '(': case ',': if ( tr.charAt(++x) != '(' ) totip[x] = --s; --x; break; } totip[x] = s; }
+
+ are = new ArrayList<Double>(n); std = new BitSet(n); u = 0;
while ( (u=tr.indexOf(")", ++u)) != last ) {
//# second tr rooting: ((ST1),(ST2),(ST3)); => ((ST1),((ST2),(ST3))); ##########################################################
//# | |
//# u last
if ( u == sup ) {
last = tro.lastIndexOf(")"); tr = tro.insert(last, ')'); v = apc(tr, apc(tr, last)); tr = tr.insert(++v, '(');
- ++last; //# NOTE: closing parenthesis at index 'last' should not be considered for REQ calculations
- /*if ( verbose ) System.out.println("# " + tr.toString());*/
+ ++last; //# NOTE: closing parenthesis at index 'last' should not be considered for REQ calculations
+ x = tr.length(); s = (short) n; while ( --x >= 0 ) { switch ( tr.charAt(x) ) { case '(': case ',': if ( tr.charAt(++x) != '(' ) totip[x] = --s; --x; break; } totip[x] = s; }
sup = 0; --u; continue;
}
//# parsing every internal branch e at index u in order to obtain lbl(STa) lbl(STb) | lbl(STc) lbl(T)-lbl(STa U STb U STc) ########
+ //# v w 'u'
+ //# | | |
//# ( ... (((STa),(STb)),STc) ... );
- //# || | | |
- //# xv w 'u' y
- v = aop(tr, u); w = apc(tr, u); x = cop(tr, u); y = ccp(tr, u); stdts = new TreeSet<Short>(ts);
- sta = new short[j=(split=tr.substring(v+1, w).replaceAll("[,\\(\\)]+", " ").trim().split(" ")).length]; for (String s: split) stdts.remove(Short.valueOf(sta[--j]=Short.parseShort(s)));
- stb = new short[j=(split=tr.substring(++w, u).replaceAll("[,\\(\\)]+", " ").trim().split(" ")).length]; for (String s: split) stdts.remove(Short.valueOf(stb[--j]=Short.parseShort(s)));
- stc = new short[j=(split=((u+1 != y) ? tr.substring(u+2, y) : tr.substring(++x, --v)).replaceAll("[,\\(\\)]+", " ").trim().split(" ")).length]; for (String s: split) stdts.remove(Short.valueOf(stc[--j]=Short.parseShort(s)));
- std = new short[j=stdts.size()]; for (Short si: stdts) std[--j] = si.shortValue();
- //# comparing every quartet ab|cd with the topology induced by dm, where a,b,c,d belongs to STa,STb,STc,STd, respectively ##########
- up = 0; //Arrays.sort(sta); Arrays.sort(stb); Arrays.sort(stb);
- for (short a: sta) {
- for (short b: stb) { ab = (a>b)?dm[a][b]:dm[b][a];
- for (short c: stc) { ac = (a>c)?dm[a][c]:dm[c][a]; bc = (b>c)?dm[b][c]:dm[c][b];
- for (short d: std) up += ( (abcd=ab+((c>d)?dm[c][d]:dm[d][c])) < ac+((b>d)?dm[b][d]:dm[d][b]) && abcd < bc+((a>d)?dm[a][d]:dm[d][a]) ) ? 1 : 0;
+ //# | |
+ //# x y
+ //# | |
+ //# ( ... (STc,((STa),(STb))) ... );
+ //# | | |
+ //# v w 'u'
+ v = aop(tr, u); w = apc(tr, u); x = cop(tr, u); y = ccp(tr, u);
+ as = totip[++v]; --v; ae = totip[w]; //# NOTE: sta = tax list inside tr.substring(v+1, w)
+ bs = totip[++w]; be = totip[u]; //# NOTE: stb = tax list inside tr.substring(++w, u)
+ if ( ++u != y ) { cs = totip[++u]; --u; ce = totip[y]; } else { cs = totip[++x]; ce = totip[--v]; } //# NOTE: stc = tax list inside ( u+1 != y ) ? tr.substring(u+2, y) : tr.substring(++x, --v)
+ --u;
+ //if ( ae-as > be-bs ) { s = ae; ae = be; be = s; s = as; as = bs; bs = s; }
+ std.set(0, n);
+ sta = new short[(i=ae-as)]; while ( --ae >= as ) std.clear((sta[--i]=tip[ae]));
+ stb = new short[(i=be-bs)]; while ( --be >= bs ) std.clear((stb[--i]=tip[be]));
+ stc = new short[(i=ce-cs)]; while ( --ce >= cs ) std.clear((stc[--i]=tip[ce]));
+ //# comparing every quartet ab|cd with the topology induced by dm, where a,b,c,d belong to STa,STb,STc,STd, respectively ###########
+ up = 0; d = -1;
+ while ( (d=std.nextSetBit(++d)) >= 0 ) {
+ for (short a: sta) { ad = (a>d)?dm[a][d]:dm[d][a];
+ for (short b: stb) { ab = (a>b)?dm[a][b]:dm[b][a]; bd = (b>d)?dm[b][d]:dm[d][b];
+ for (short c: stc) up += ( (abcd=ab+((c>d)?dm[c][d]:dm[d][c])) < ((a>c)?dm[a][c]:dm[c][a])+bd && abcd < ad+((b>c)?dm[b][c]:dm[c][b]) ) ? 0b1 : 0;
}
}
}
//# storing re value inside are #####################################################################################################
- are.add(Double.valueOf(re=up/(dn=((double)sta.length)*((double)stb.length)*((double)stc.length)*((double)std.length))));
- if ( verbose ) System.out.println((Arrays.toString(sta) + Arrays.toString(stb) + Arrays.toString(stc) + Arrays.toString(std)).replaceAll(" ","") + " Re=" + df.format(re) + " (" + ((long)up) + "/" + ((long)dn) + ")");
+ are.add(Double.valueOf(re=up/(dn=((double)(sta.length))*((double)(stb.length))*((double)(stc.length))*((double)std.cardinality()))));
+ if ( verbose ) {
+ Arrays.sort(sta); Arrays.sort(stb); Arrays.sort(stc);
+ System.out.println((Arrays.toString(sta) + Arrays.toString(stb) + Arrays.toString(stc) + std.toString().replace('{','[').replace('}',']')).replaceAll(" ","")
+ + String.format(" Re=%.3f", re) + " (" + ((long)up) + "/" + ((long)dn) + ")");
+ }
}
//### writing re values into nwk #####################################################################################################################################################
- u = -1; v = -1; while ( (u=nwk.indexOf(")", ++u)) != -1 ) { if ( nwk.charAt(u+1) == ';' ) break; nwk = nwk.insert(u+1, df.format(are.get(++v).doubleValue())); }
- out = new BufferedWriter(new FileWriter(new File(args[2]))); out.write(nwk.toString()); out.close(); //System.out.println(nwk.toString());
+ u = 0; v = -1; while ( (u=nwk.indexOf(")", u)) != -1 ) if ( nwk.charAt(++u) != ';' ) nwk = nwk.insert(u, String.format(Locale.US, "%.3f", are.get(++v)));
+ out = new BufferedWriter(new FileWriter(new File(args[2]))); out.write(nwk.toString()); out.newLine(); out.close();
if ( verbose ) System.out.println("output tree written into " + args[2]);
}
- //## containing opening parenthesis (cop): returns x from t and i, where i is a closing parenthesis index
+ //## containing opening parenthesis (cop): returns x from t and i, where i is a closing parenthesis index ################
// t = ( ... ((STa),(STb)) ... );
// | |
// x i
- final static int cop (StringBuilder t, int i) { short p = 0; ++p; int x = i; while ( --x >= 0 ) switch ( t.charAt(x) ) { case ')': ++p; continue; case '(': if ( --p < 0 ) return x; } return -1; }
+ final static int cop (final StringBuilder t, final int i) { short p = 0b1; int x = i; while ( --x >= 0 ) switch ( t.charAt(x) ) { case ')': ++p; continue; case '(': if ( --p < 0 ) return x; } return -1; }
- //## containing closing parenthesis (ccp): returns x from t and i, where i is a closing parenthesis index
+ //## containing closing parenthesis (ccp): returns x from t and i, where i is a closing parenthesis index ################
// t = ( ... ((STa),(STb)) ... );
// | |
// i x
- final static int ccp (StringBuilder t, int i) { short p = 0; int x = i; while ( true ) switch ( t.charAt(++x) ) { case ';': return -1; case '(': ++p; continue; case ')': if ( --p < 0 ) return x; } }
+ final static int ccp (final StringBuilder t, final int i) { short p = 0; int x = i; while ( true ) switch ( t.charAt(++x) ) { case ';': return -1; case '(': ++p; continue; case ')': if ( --p < 0 ) return x; } }
- //## associated opening parenthesis (aop): returns x from t and i, where i is a closing parenthesis index
+ //## associated opening parenthesis (aop): returns x from t and i, where i is a closing parenthesis index ################
// t = ( ... ((STa),(STb)) ... );
// | |
// x i
- final static int aop (StringBuilder t, int i) { short p = 0; ++p; int x = i; while ( --x >= 0 ) switch ( t.charAt(x) ) { case ')': ++p; continue; case '(': if ( --p == 0 ) return x; } return -1; }
+ final static int aop (final StringBuilder t, final int i) { short p = 0b1; int x = i; while ( --x >= 0 ) switch ( t.charAt(x) ) { case ')': ++p; continue; case '(': if ( --p == 0 ) return x; } return -1; }
- //## associated partitionning comma (apc): returns x from t and i, where i is a closing parenthesis index
+ //## associated partitionning comma (apc): returns x from t and i, where i is a closing parenthesis index ################
// t = ( ... ((STa),(STb)) ... );
// | |
// x i
- final static int apc (StringBuilder t, int i) { short p = 0; int x = i; while ( --x >= 0 ) switch ( t.charAt(x) ) { case ')': ++p; continue; case '(': --p; continue; case ',': if ( p == 0 ) return x; } return -1; }
+ final static int apc (final StringBuilder t, final int i) { short p = 0; int x = i; while ( --x >= 0 ) switch ( t.charAt(x) ) { case ')': ++p; continue; case '(': --p; continue; case ',': if ( p == 0 ) return x; } return -1; }
}
+