//*****************************************************************************
//
// File: MaximumParsimonyBnbSmp.java
// Package: edu.rit.compbio.phyl
// Unit: Class edu.rit.compbio.phyl.MaximumParsimonyBnbSmp
//
// This Java source file is copyright (C) 2008 by Alan Kaminsky. All rights
// reserved. For further information, contact the author, Alan Kaminsky, at
// ark@cs.rit.edu.
//
// This Java source file is part of the Parallel Java Library ("PJ"). PJ 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.
//
// PJ 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.
//
// Linking this library statically or dynamically with other modules is making a
// combined work based on this library. Thus, the terms and conditions of the
// GNU General Public License cover the whole combination.
//
// As a special exception, the copyright holders of this library give you
// permission to link this library with independent modules to produce an
// executable, regardless of the license terms of these independent modules, and
// to copy and distribute the resulting executable under terms of your choice,
// provided that you also meet, for each linked independent module, the terms
// and conditions of the license of that module. An independent module is a
// module which is not derived from or based on this library. If you modify this
// library, you may extend this exception to your version of the library, but
// you are not obligated to do so. If you do not wish to do so, delete this
// exception statement from your version.
//
// A copy of the GNU General Public License is provided in the file gpl.txt. You
// may also obtain a copy of the GNU General Public License on the World Wide
// Web at http://www.gnu.org/licenses/gpl.html.
//
//******************************************************************************
package edu.rit.compbio.phyl;
import edu.rit.pj.reduction.IntegerOp;
import edu.rit.pj.reduction.SharedInteger;
/**
* Class MaximumParsimonyBnbSmp provides an SMP parallel algorithm for maximum
* parsimony phylogenetic tree construction using branch-and-bound search.
*
* Class MaximumParsimonyBnbSmp is designed to be used in an SMP parallel
* program that runs in one process with multiple threads. Each thread has its
* own MaximumParsimonyBnbSmp instance. Each thread uses its own
* MaximumParsimonyBnbSmp instance to search different sections of the search
* graph concurrently.
*
* The process has a shared variable, bound, that holds the best
* parsimony score found so far (i.e., the bound for branch-and-bound search).
* The bound variable is an instance of class {@linkplain
* edu.rit.pj.reduction.SharedInteger}. All the MaximumParsimonyBnbSmp instances
* within the process share the same bound variable. Whenever one thread
* finds a phylogenetic tree with a better parsimony score, it notifies all the
* threads by updating the bound variable.
*
* To perform a search, the process must:
*
* -
* Call the static createBoundVariable() method to create the shared
* bound variable.
*
* Then each thread in the process must:
*
* -
* Create its own instance of class MaximumParsimonyBnbSmp, passing in a
* {@linkplain DnaSequenceList} of the DNA sequences in the tree, a reference to
* the process's shared bound variable, and a {@linkplain
* MaximumParsimonyResults} object to hold the search results.
*
-
* Call the findTrees() method one or more times, each time indicating
* a section of the search graph to search. The results of searching that
* section are accumulated into the {@linkplain MaximumParsimonyResults} object
* specified to the constructor.
*
*
* Note: Class MaximumParsimonyBnbSmp is not multiple thread safe; it is
* intended to be used as a per-thread variable.
*
* @author Alan Kaminsky
* @version 21-Nov-2008
*/
public class MaximumParsimonyBnbSmp
{
// Exported static operations.
/**
* Create a new shared bound variable.
*
* @param initialBound Initial bound for branch-and-bound search.
*/
public static SharedInteger createBoundVariable
(int initialBound)
{
return new SharedInteger (initialBound);
}
// Hidden data members.
// List of DNA sequences with which to construct trees.
private DnaSequenceList seqList;
// Shared bound variable.
private SharedInteger bound;
// For holding search results.
private MaximumParsimonyResults results;
// Length of each DNA sequence.
private int L;
// Number of DNA sequences.
private int N;
// Tree capacity.
private int C;
// Number of absent states as each DNA sequence is added.
private int[] absentStates;
// Stack of DNA sequence trees.
private DnaSequenceTree[] treeStack;
// Stack of auxiliary DNA sequence arrays.
DnaSequence[][] seqArrayStack;
// Tree signature currently being searched.
private int[] signature;
// Extra padding to avert cache interference.
private long p0, p1, p2, p3, p4, p5, p6, p7;
private long p8, p9, pa, pb, pc, pd, pe, pf;
// Exported constructors.
/**
* Construct a new maximum parsimony phylogenetic tree construction
* algorithm object.
*
* @param seqList DNA sequence list.
* @param bound Shared bound variable.
* @param results Object in which to store the results.
*/
public MaximumParsimonyBnbSmp
(DnaSequenceList seqList,
SharedInteger bound,
MaximumParsimonyResults results)
{
// Record parameters.
this.seqList = seqList;
this.bound = bound;
this.results = results;
// Initialize.
L = seqList.seq(0).length();
N = seqList.length();
C = 2*N - 1;
// Compute number of absent states as each DNA sequence is added.
absentStates = seqList.countAbsentStates();
// Set up stack of DNA sequence trees.
treeStack = new DnaSequenceTree [N];
for (int i = 0; i < N; ++ i)
{
treeStack[i] = new DnaSequenceTree (C);
}
// Initialize DNA sequence tree at first level of the search graph.
treeStack[0].add (0, seqList.seq(0));
// Set up stack of auxiliary DNA sequence arrays.
seqArrayStack = new DnaSequence [N] [];
for (int i = 0; i < N; ++ i)
{
DnaSequence[] seqArray = new DnaSequence [i];
seqArrayStack[i] = seqArray;
for (int j = 0; j < i; ++ j)
{
seqArray[j] = new DnaSequence (L);
}
}
// Set up tree signature.
signature = new int [N+32]; // Extra padding
}
// Exported operations.
/**
* Find the maximum parsimony phylogenetic tree(s) in the given section of
* the search graph. The DNA sequence list was specified to the constructor.
*
* The search will commence at level L of the search graph, 0 ≤
* L ≤ N−1, where N is the number of sequences
* in the DNA sequence list. Of the (2L−1)!! vertices at level
* L, the search will start at the V1-th such
* vertex and end at the V2-th such vertex, 0 ≤
* V1 ≤ V2 ≤
* (2L−1)!! − 1. All search graph vertices at
* and below vertices V1 through V2
* inclusive will be searched.
*
* The results are accumulated into the {@linkplain MaximumParsimonyResults}
* object specified to the constructor. The findTrees() method will
* only find trees whose parsimony scores are less than or equal to the
* value of the bound variable specified to the constructor or the
* best bound found thereafter, whichever is smaller.
*
* @param startLevel
* L, the level of the search graph at which to commence the
* search.
* @param vertex1
* V1, the search graph vertex at level L at
* which to start the search.
* @param vertex2
* V2, the search graph vertex at level L at
* which to end the search.
*/
public void findTrees
(int startLevel,
int vertex1,
int vertex2)
{
// Initialize tree signature as specified by and .
signature[0] = 0;
int q = vertex1;
for (int i = startLevel; i > 0; -- i)
{
int d = 2*i - 1;
signature[i] = q % d - 1;
q = q / d;
}
for (int i = startLevel + 1; i < N; ++ i)
{
signature[i] = -1;
}
// Traverse remaining levels of the search graph.
int level = 1;
boolean done = false;
results.reduceScore (bound.get());
while (! done)
{
DnaSequenceTree prevTree = treeStack[level-1];
// If we have reached the bottom of the search graph, we have a
// tentative solution.
if (level == N)
{
int tentativeScore = prevTree.seq (prevTree.root()) .score();
// Update best solution's score to reflect tentative solution's
// score.
int updatedScore =
bound.reduce (tentativeScore, IntegerOp.MINIMUM);
// If best solution's score is better than that of previous
// solutions, discard previous solutions.
results.reduceScore (updatedScore);
// Record tentative solution.
results.add (signature, tentativeScore);
// Go to previous level.
-- level;
if (level == startLevel)
{
++ vertex1;
if (vertex1 > vertex2) done = true;
}
}
// If there are no more positions to try at this level, reset
// position at this level and go to previous level.
else if (signature[level] == 2*(level - 1))
{
signature[level] = -1;
-- level;
if (level == startLevel)
{
++ vertex1;
if (vertex1 > vertex2) done = true;
}
}
// If there are more positions to try at this level, add the DNA
// sequence to the tree at the next position and do
// branch-and-bound.
else
{
++ signature[level];
DnaSequenceTree currTree = treeStack[level];
currTree.copy (prevTree);
int tip = currTree.add (signature[level], seqList.seq(level));
int partialScore =
FitchParsimony.updateScore
(currTree, tip, seqArrayStack[level]);
// If partial parsimony score plus number of absent states in
// the remaining levels is less than or equal to the best
// solution's score, go to the next level (branch), otherwise
// try the next choice at this level (bound).
if
((level <= startLevel) ||
(partialScore + absentStates[level] <= bound.get()))
{
++ level;
}
}
}
}
}