gsdi work
[jalview.git] / forester / java / src / org / forester / sdi / GSDI.java
index 6947910..a10a085 100644 (file)
@@ -27,10 +27,10 @@ package org.forester.sdi;
 
 import java.util.ArrayList;
 import java.util.HashMap;
+import java.util.HashSet;
 import java.util.List;
 import java.util.Map;
-import java.util.SortedSet;
-import java.util.TreeSet;
+import java.util.Set;
 
 import org.forester.phylogeny.Phylogeny;
 import org.forester.phylogeny.PhylogenyNode;
@@ -63,43 +63,15 @@ import org.forester.util.ForesterUtil;
  */
 public final class GSDI extends SDI {
 
-    private final HashMap<PhylogenyNode, Integer> _transversal_counts;
-    private final boolean                         _most_parsimonious_duplication_model;
-    private final boolean                         _strip_gene_tree;
-    private final boolean                         _strip_species_tree;
-    private int                                   _speciation_or_duplication_events_sum;
-    private int                                   _speciations_sum;
-    private final List<PhylogenyNode>             _stripped_gene_tree_nodes;
-    private final List<PhylogenyNode>             _stripped_species_tree_nodes;
-    private final SortedSet<PhylogenyNode>        _mapped_species_tree_nodes;
+    private final boolean             _most_parsimonious_duplication_model;
+    private final boolean             _strip_gene_tree;
+    private final boolean             _strip_species_tree;
+    private int                       _speciation_or_duplication_events_sum;
+    private int                       _speciations_sum;
+    private final List<PhylogenyNode> _stripped_gene_tree_nodes;
+    private final List<PhylogenyNode> _stripped_species_tree_nodes;
+    private final Set<PhylogenyNode>  _mapped_species_tree_nodes;
 
-    /**
-     * Constructor which sets the gene tree and the species tree to be compared.
-     * species_tree is the species tree to which the gene tree gene_tree will be
-     * compared to - with method "infer(boolean)". Both Trees must be completely
-     * binary and rooted. The actual inference is accomplished with method
-     * "infer(boolean)". The mapping cost L can then be calculated with method
-     * "computeMappingCost()".
-     * <p>
-     * 
-     * @see #infer(boolean)
-     * @see SDI#computeMappingCostL()
-     * @param gene_tree
-     *            reference to a rooted gene tree to which assign duplication vs
-     *            speciation, must have species names in the species name fields
-     *            for all external nodes
-     * @param species_tree
-     *            reference to a rooted binary species tree which might get
-     *            stripped in the process, must have species names in the
-     *            species name fields for all external nodes
-     * 
-     * @param most_parsimonious_duplication_model
-     *            set to true to assign nodes as speciations which would
-     *            otherwise be assiged as unknown because of polytomies in the
-     *            species tree.
-     * @throws SdiException 
-     * 
-     */
     public GSDI( final Phylogeny gene_tree,
                  final Phylogeny species_tree,
                  final boolean most_parsimonious_duplication_model,
@@ -109,16 +81,15 @@ public final class GSDI extends SDI {
         _speciation_or_duplication_events_sum = 0;
         _speciations_sum = 0;
         _most_parsimonious_duplication_model = most_parsimonious_duplication_model;
-        _transversal_counts = new HashMap<PhylogenyNode, Integer>();
         _duplications_sum = 0;
         _strip_gene_tree = strip_gene_tree;
         _strip_species_tree = strip_species_tree;
         _stripped_gene_tree_nodes = new ArrayList<PhylogenyNode>();
         _stripped_species_tree_nodes = new ArrayList<PhylogenyNode>();
-        _mapped_species_tree_nodes = new TreeSet<PhylogenyNode>();
+        _mapped_species_tree_nodes = new HashSet<PhylogenyNode>();
         getSpeciesTree().preOrderReId();
         linkNodesOfG();
-        geneTreePostOrderTraversal( getGeneTree().getRoot() );
+        geneTreePostOrderTraversal();
     }
 
     GSDI( final Phylogeny gene_tree, final Phylogeny species_tree, final boolean most_parsimonious_duplication_model )
@@ -126,78 +97,58 @@ public final class GSDI extends SDI {
         this( gene_tree, species_tree, most_parsimonious_duplication_model, false, false );
     }
 
-    private final Event createDuplicationEvent() {
-        final Event event = Event.createSingleDuplicationEvent();
-        ++_duplications_sum;
-        return event;
-    }
-
-    private final Event createSingleSpeciationOrDuplicationEvent() {
-        final Event event = Event.createSingleSpeciationOrDuplicationEvent();
-        ++_speciation_or_duplication_events_sum;
-        return event;
-    }
-
-    private final Event createSpeciationEvent() {
-        final Event event = Event.createSingleSpeciationEvent();
-        ++_speciations_sum;
-        return event;
-    }
-
     // s is the node on the species tree g maps to.
     private final void determineEvent( final PhylogenyNode s, final PhylogenyNode g ) {
-        Event event = null;
-        // Determine how many children map to same node as parent.
-        int sum_g_childs_mapping_to_s = 0;
-        for( final PhylogenyNodeIterator iter = g.iterateChildNodesForward(); iter.hasNext(); ) {
-            if ( iter.next().getLink() == s ) {
-                ++sum_g_childs_mapping_to_s;
-            }
+        boolean oyako = false;
+        if ( ( g.getChildNode1().getLink() == s ) || ( g.getChildNode2().getLink() == s ) ) {
+            oyako = true;
         }
-        // Determine the sum of traversals.
-        int traversals_sum = 0;
-        int max_traversals = 0;
-        PhylogenyNode max_traversals_node = null;
-        if ( !s.isExternal() ) {
-            for( final PhylogenyNodeIterator iter = s.iterateChildNodesForward(); iter.hasNext(); ) {
-                final PhylogenyNode current_node = iter.next();
-                final int traversals = getTraversalCount( current_node );
-                traversals_sum += traversals;
-                if ( traversals > max_traversals ) {
-                    max_traversals = traversals;
-                    max_traversals_node = current_node;
-                }
+        if ( g.getLink().getNumberOfDescendants() == 2 ) {
+            if ( oyako ) {
+                g.getNodeData().setEvent( createDuplicationEvent() );
             }
-        }
-        // System.out.println( " sum=" + traversals_sum );
-        // System.out.println( " max=" + max_traversals );
-        // System.out.println( " m=" + sum_g_childs_mapping_to_s );
-        if ( sum_g_childs_mapping_to_s > 0 ) {
-            if ( traversals_sum == 2 ) {
-                event = createDuplicationEvent();
+            else {
+                g.getNodeData().setEvent( createSpeciationEvent() );
             }
-            else if ( traversals_sum > 2 ) {
-                if ( max_traversals <= 1 ) {
-                    if ( _most_parsimonious_duplication_model ) {
-                        event = createSpeciationEvent();
+        }
+        else {
+            if ( oyako ) {
+                final Set<PhylogenyNode> set = new HashSet<PhylogenyNode>();
+                for( PhylogenyNode n : g.getChildNode1().getAllExternalDescendants() ) {
+                    n = n.getLink();
+                    while ( n.getParent() != s ) {
+                        n = n.getParent();
+                        if ( n.isRoot() ) {
+                            break;
+                        }
                     }
-                    else {
-                        event = createSingleSpeciationOrDuplicationEvent();
+                    set.add( n );
+                }
+                boolean multiple = false;
+                for( PhylogenyNode n : g.getChildNode2().getAllExternalDescendants() ) {
+                    n = n.getLink();
+                    while ( n.getParent() != s ) {
+                        n = n.getParent();
+                        if ( n.isRoot() ) {
+                            break;
+                        }
                     }
+                    if ( set.contains( n ) ) {
+                        multiple = true;
+                        break;
+                    }
+                }
+                if ( multiple ) {
+                    g.getNodeData().setEvent( createDuplicationEvent() );
                 }
                 else {
-                    event = createDuplicationEvent();
-                    _transversal_counts.put( max_traversals_node, 1 );
+                    g.getNodeData().setEvent( createSingleSpeciationOrDuplicationEvent() );
                 }
             }
             else {
-                event = createDuplicationEvent();
+                g.getNodeData().setEvent( createSpeciationEvent() );
             }
         }
-        else {
-            event = createSpeciationEvent();
-        }
-        g.getNodeData().setEvent( event );
     }
 
     /**
@@ -209,65 +160,51 @@ public final class GSDI extends SDI {
      * the species tree must be labeled in preorder.
      * <p>
      * 
-     * @param g
-     *            starting node of a gene tree - normally the root
      */
-    final void geneTreePostOrderTraversal( final PhylogenyNode g ) {
-        if ( !g.isExternal() ) {
-            for( final PhylogenyNodeIterator iter = g.iterateChildNodesForward(); iter.hasNext(); ) {
-                geneTreePostOrderTraversal( iter.next() );
-            }
-            final PhylogenyNode[] linked_nodes = new PhylogenyNode[ g.getNumberOfDescendants() ];
-            for( int i = 0; i < linked_nodes.length; ++i ) {
-                if ( g.getChildNode( i ).getLink() == null ) {
-                    System.out.println( "link is null for " + g.getChildNode( i ) );
-                    System.exit( -1 );
+    final void geneTreePostOrderTraversal() {
+        for( final PhylogenyNodeIterator it = getGeneTree().iteratorPostorder(); it.hasNext(); ) {
+            final PhylogenyNode g = it.next();
+            if ( g.isInternal() ) {
+                PhylogenyNode s1 = g.getChildNode1().getLink();
+                PhylogenyNode s2 = g.getChildNode2().getLink();
+                while ( s1 != s2 ) {
+                    if ( s1.getId() > s2.getId() ) {
+                        s1 = s1.getParent();
+                    }
+                    else {
+                        s2 = s2.getParent();
+                    }
                 }
-                linked_nodes[ i ] = g.getChildNode( i ).getLink();
-            }
-            final int[] min_max = obtainMinMaxIdIndices( linked_nodes );
-            int min_i = min_max[ 0 ];
-            int max_i = min_max[ 1 ];
-            // initTransversalCounts();
-            while ( linked_nodes[ min_i ] != linked_nodes[ max_i ] ) {
-                increaseTraversalCount( linked_nodes[ max_i ] );
-                linked_nodes[ max_i ] = linked_nodes[ max_i ].getParent();
-                final int[] min_max_ = obtainMinMaxIdIndices( linked_nodes );
-                min_i = min_max_[ 0 ];
-                max_i = min_max_[ 1 ];
+                g.setLink( s1 );
+                determineEvent( s1, g );
             }
-            final PhylogenyNode s = linked_nodes[ max_i ];
-            g.setLink( s );
-            // Determines whether dup. or spec.
-            determineEvent( s, g );
-            // _transversal_counts.clear();
         }
     }
 
-    public final int getSpeciationOrDuplicationEventsSum() {
-        return _speciation_or_duplication_events_sum;
+    private final Event createDuplicationEvent() {
+        final Event event = Event.createSingleDuplicationEvent();
+        ++_duplications_sum;
+        return event;
     }
 
-    public final int getSpeciationsSum() {
-        return _speciations_sum;
+    private final Event createSingleSpeciationOrDuplicationEvent() {
+        final Event event = Event.createSingleSpeciationOrDuplicationEvent();
+        ++_speciation_or_duplication_events_sum;
+        return event;
     }
 
-    private final int getTraversalCount( final PhylogenyNode node ) {
-        if ( _transversal_counts.containsKey( node ) ) {
-            return _transversal_counts.get( node );
-        }
-        return 0;
+    private final Event createSpeciationEvent() {
+        final Event event = Event.createSingleSpeciationEvent();
+        ++_speciations_sum;
+        return event;
     }
 
-    private final void increaseTraversalCount( final PhylogenyNode node ) {
-        if ( _transversal_counts.containsKey( node ) ) {
-            _transversal_counts.put( node, _transversal_counts.get( node ) + 1 );
-        }
-        else {
-            _transversal_counts.put( node, 1 );
-        }
-        // System.out.println( "count for node " + node.getID() + " is now "
-        // + getTraversalCount( node ) );
+    public final int getSpeciationOrDuplicationEventsSum() {
+        return _speciation_or_duplication_events_sum;
+    }
+
+    public final int getSpeciationsSum() {
+        return _speciations_sum;
     }
 
     /**
@@ -277,34 +214,11 @@ public final class GSDI extends SDI {
      * 
      */
     @Override
-    //    final void linkNodesOfG() {
-    //        final HashMap<Taxonomy, PhylogenyNode> speciestree_ext_nodes = createTaxonomyToNodeMap();
-    //        if ( _strip_gene_tree ) {
-    //            stripGeneTree( speciestree_ext_nodes );
-    //            if ( ( _gene_tree == null ) || ( _gene_tree.getNumberOfExternalNodes() < 2 ) ) {
-    //                throw new IllegalArgumentException( "species tree does not contain any"
-    //                        + " nodes matching species in the gene tree" );
-    //            }
-    //        }
-    //        // Retrieve the reference to the PhylogenyNode with a matching species.
-    //        for( final PhylogenyNodeIterator iter = _gene_tree.iteratorExternalForward(); iter.hasNext(); ) {
-    //            final PhylogenyNode g = iter.next();
-    //            if ( !g.getNodeData().isHasTaxonomy() ) {
-    //                throw new IllegalArgumentException( "gene tree node " + g + " has no taxonomic data" );
-    //            }
-    //            final PhylogenyNode s = speciestree_ext_nodes.get( g.getNodeData().getTaxonomy() );
-    //            if ( s == null ) {
-    //                throw new IllegalArgumentException( "species " + g.getNodeData().getTaxonomy()
-    //                        + " not present in species tree" );
-    //            }
-    //            g.setLink( s );
-    //        }
-    //    }
     final void linkNodesOfG() throws SdiException {
         final Map<String, PhylogenyNode> species_to_node_map = new HashMap<String, PhylogenyNode>();
         final List<PhylogenyNode> species_tree_ext_nodes = new ArrayList<PhylogenyNode>();
         final TaxonomyComparisonBase tax_comp_base = determineTaxonomyComparisonBase( _gene_tree );
-        System.out.println( "comp base is: " + tax_comp_base );
+        // System.out.println( "comp base is: " + tax_comp_base );
         // Stringyfied taxonomy is the key, node is the value.
         for( final PhylogenyNodeIterator iter = _species_tree.iteratorExternalForward(); iter.hasNext(); ) {
             final PhylogenyNode s = iter.next();
@@ -352,7 +266,7 @@ public final class GSDI extends SDI {
                     else {
                         g.setLink( s );
                         _mapped_species_tree_nodes.add( s );
-                        System.out.println( "setting link of " + g + " to " + s );
+                        //  System.out.println( "setting link of " + g + " to " + s );
                     }
                 }
             }
@@ -363,75 +277,18 @@ public final class GSDI extends SDI {
             }
         }
         if ( _strip_species_tree ) {
-            for( PhylogenyNode x : _mapped_species_tree_nodes ) {
-                System.out.println( ">>" + x );
-            }
             for( final PhylogenyNode s : species_tree_ext_nodes ) {
-                System.out.print( ">>>>>>>>>" + s );
                 if ( !_mapped_species_tree_nodes.contains( s ) ) {
                     _species_tree.deleteSubtree( s, true );
-                    System.out.println( " DELETING" );
                 }
-                else {
-                    System.out.println();
-                }
-            }
-            for( PhylogenyNode x : _mapped_species_tree_nodes ) {
-                System.out.println( ">>" + x );
             }
         }
     }
 
-    public SortedSet<PhylogenyNode> getMappedExternalSpeciesTreeNodes() {
+    public Set<PhylogenyNode> getMappedExternalSpeciesTreeNodes() {
         return _mapped_species_tree_nodes;
     }
 
-    //    final private HashMap<Taxonomy, PhylogenyNode> createTaxonomyToNodeMap() {
-    //        final HashMap<Taxonomy, PhylogenyNode> speciestree_ext_nodes = new HashMap<Taxonomy, PhylogenyNode>();
-    //        for( final PhylogenyNodeIterator iter = _species_tree.iteratorLevelOrder(); iter.hasNext(); ) {
-    //            final PhylogenyNode n = iter.next();
-    //            if ( n.getNodeData().isHasTaxonomy() ) {
-    //                if ( speciestree_ext_nodes.containsKey( n.getNodeData().getTaxonomy() ) ) {
-    //                    throw new IllegalArgumentException( "taxonomy [" + n.getNodeData().getTaxonomy()
-    //                            + "] is not unique in species phylogeny" );
-    //                }
-    //                speciestree_ext_nodes.put( n.getNodeData().getTaxonomy(), n );
-    //            }
-    //        }
-    //        return speciestree_ext_nodes;
-    //    }
-    //    private final void stripGeneTree( final HashMap<Taxonomy, PhylogenyNode> speciestree_ext_nodes ) {
-    //        //  final Set<PhylogenyNode> to_delete = new HashSet<PhylogenyNode>();
-    //        for( final PhylogenyNodeIterator iter = _gene_tree.iteratorExternalForward(); iter.hasNext(); ) {
-    //            final PhylogenyNode g = iter.next();
-    //            if ( !g.getNodeData().isHasTaxonomy() ) {
-    //                throw new IllegalArgumentException( "gene tree node " + g + " has no taxonomic data" );
-    //            }
-    //            if ( !speciestree_ext_nodes.containsKey( g.getNodeData().getTaxonomy() ) ) {
-    //                _stripped_gene_tree_nodes.add( g );
-    //            }
-    //        }
-    //        for( final PhylogenyNode n : _stripped_gene_tree_nodes ) {
-    //            _gene_tree.deleteSubtree( n, true );
-    //        }
-    //    }
-    //    private final void stripGeneTree2( final HashMap<Taxonomy, PhylogenyNode> speciestree_ext_nodes ) {
-    //        //  final Set<PhylogenyNode> to_delete = new HashSet<PhylogenyNode>();
-    //        for( final PhylogenyNodeIterator iter = _gene_tree.iteratorExternalForward(); iter.hasNext(); ) {
-    //            final PhylogenyNode g = iter.next();
-    //            if ( !g.getNodeData().isHasTaxonomy() ) {
-    //                _stripped_gene_tree_nodes.add( g );
-    //            }
-    //            else {
-    //                if ( !speciestree_ext_nodes.containsKey( g.getNodeData().getTaxonomy() ) ) {
-    //                    _stripped_gene_tree_nodes.add( g );
-    //                }
-    //            }
-    //        }
-    //        for( final PhylogenyNode n : _stripped_gene_tree_nodes ) {
-    //            _gene_tree.deleteSubtree( n, true );
-    //        }
-    //    }
     public static TaxonomyComparisonBase determineTaxonomyComparisonBase( final Phylogeny gene_tree ) {
         int with_id_count = 0;
         int with_code_count = 0;
@@ -495,96 +352,4 @@ public final class GSDI extends SDI {
         sb.append( "mapping cost L                     : " + computeMappingCostL() );
         return sb.toString();
     }
-
-    static final int[] obtainMinMaxIdIndices( final PhylogenyNode[] linked_nodes ) {
-        int max_i = 0;
-        int min_i = 0;
-        int max_i_id = -Integer.MAX_VALUE;
-        int min_i_id = Integer.MAX_VALUE;
-        for( int i = 0; i < linked_nodes.length; ++i ) {
-            final int id_i = linked_nodes[ i ].getId();
-            if ( id_i > max_i_id ) {
-                max_i = i;
-                max_i_id = linked_nodes[ max_i ].getId();
-            }
-            if ( id_i < min_i_id ) {
-                min_i = i;
-                min_i_id = linked_nodes[ min_i ].getId();
-            }
-        }
-        return new int[] { min_i, max_i };
-    }
-    /**
-     * Updates the mapping function M after the root of the gene tree has been
-     * moved by one branch. It calculates M for the root of the gene tree and
-     * one of its two children.
-     * <p>
-     * To be used ONLY by method "SDIunrooted.fastInfer(Phylogeny,Phylogeny)".
-     * <p>
-     * (Last modfied: )
-     * 
-     * @param prev_root_was_dup
-     *            true if the previous root was a duplication, false otherwise
-     * @param prev_root_c1
-     *            child 1 of the previous root
-     * @param prev_root_c2
-     *            child 2 of the previous root
-     * @return number of duplications which have been assigned in gene tree
-     */
-    // int updateM( final boolean prev_root_was_dup,
-    // final PhylogenyNode prev_root_c1, final PhylogenyNode prev_root_c2 ) {
-    // final PhylogenyNode root = getGeneTree().getRoot();
-    // if ( ( root.getChildNode1() == prev_root_c1 )
-    // || ( root.getChildNode2() == prev_root_c1 ) ) {
-    // calculateMforNode( prev_root_c1 );
-    // }
-    // else {
-    // calculateMforNode( prev_root_c2 );
-    // }
-    // Event event = null;
-    // if ( prev_root_was_dup ) {
-    // event = Event.createSingleDuplicationEvent();
-    // }
-    // else {
-    // event = Event.createSingleSpeciationEvent();
-    // }
-    // root.getPhylogenyNodeData().setEvent( event );
-    // calculateMforNode( root );
-    // return getDuplications();
-    // } // updateM( boolean, PhylogenyNode, PhylogenyNode )
-    // Helper method for updateM( boolean, PhylogenyNode, PhylogenyNode )
-    // Calculates M for PhylogenyNode n, given that M for the two children
-    // of n has been calculated.
-    // (Last modified: 10/02/01)
-    // private void calculateMforNode( final PhylogenyNode n ) {
-    // if ( !n.isExternal() ) {
-    // boolean was_duplication = n.isDuplication();
-    // PhylogenyNode a = n.getChildNode1().getLink(), b = n
-    // .getChildNode2().getLink();
-    // while ( a != b ) {
-    // if ( a.getID() > b.getID() ) {
-    // a = a.getParent();
-    // }
-    // else {
-    // b = b.getParent();
-    // }
-    // }
-    // n.setLink( a );
-    // Event event = null;
-    // if ( ( a == n.getChildNode1().getLink() )
-    // || ( a == n.getChildNode2().getLink() ) ) {
-    // event = Event.createSingleDuplicationEvent();
-    // if ( !was_duplication ) {
-    // increaseDuplications();
-    // }
-    // }
-    // else {
-    // event = Event.createSingleSpeciationEvent();
-    // if ( was_duplication ) {
-    // decreaseDuplications();
-    // }
-    // }
-    // n.getPhylogenyNodeData().setEvent( event );
-    // }
-    // } // calculateMforNode( PhylogenyNode )
 }