//*****************************************************************************
//
// File:    MaximumParsimonyBnbHyb.java
// Package: edu.rit.compbio.phyl
// Unit:    Class edu.rit.compbio.phyl.MaximumParsimonyBnbHyb
//
// 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.Comm;

import edu.rit.pj.reduction.IntegerOp;

import edu.rit.pj.replica.ReplicatedInteger;

import edu.rit.util.Range;

import java.io.IOException;

/**
 * Class MaximumParsimonyBnbHyb provides a hybrid parallel algorithm for maximum
 * parsimony phylogenetic tree construction using branch-and-bound search.
 * <P>
 * Class MaximumParsimonyBnbHyb is designed to be used in a hybrid parallel
 * program that runs on multiple nodes with one process per node and multiple
 * threads per process. Each thread has its own MaximumParsimonyBnbHyb instance.
 * Each thread uses its own MaximumParsimonyBnbHyb instance to search different
 * sections of the search graph concurrently.
 * <P>
 * Each process has a shared, replicated variable, <I>bound</I>, that holds the
 * best parsimony score found so far (i.e., the bound for branch-and-bound
 * search). The <I>bound</I> variable is an instance of class {@linkplain
 * edu.rit.pj.replica.ReplicatedInteger}. All the MaximumParsimonyBnbHyb
 * instances within the process share the same <I>bound</I> variable. Whenever
 * one thread finds a phylogenetic tree with a better parsimony score, it
 * notifies all the threads and processes by updating the <I>bound</I> variable.
 * <P>
 * To perform a search, each process must:
 * <OL TYPE=1>
 * <LI>
 * Call the static <TT>createBoundVariable()</TT> method to create the shared,
 * replicated <I>bound</I> variable.
 * </OL>
 * Then each thread in each process must:
 * <OL TYPE=1 START=2>
 * <LI>
 * Create its own instance of class MaximumParsimonyBnbHyb, passing in a
 * {@linkplain DnaSequenceList} of the DNA sequences in the tree, a reference to
 * the process's shared, replicated <I>bound</I> variable, and a {@linkplain
 * MaximumParsimonyResults} object to hold the search results.
 * <P><LI>
 * Call the <TT>findTrees()</TT> 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.
 * </OL>
 * <P>
 * <I>Note:</I> Class MaximumParsimonyBnbHyb 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 MaximumParsimonyBnbHyb
	{

// Exported static operations.

	/**
	 * Create a new shared, replicated <I>bound</I> variable. The <I>bound</I>
	 * variable will use the world communicator and a message tag of
	 * <TT>Integer.MAX_VALUE</TT>. This message tag must not be used in any
	 * other messages in the program.
	 *
	 * @param  initialBound    Initial bound for branch-and-bound search.
	 */
	public static ReplicatedInteger createBoundVariable
		(int initialBound)
		{
		return createBoundVariable
			(initialBound, Comm.world(), Integer.MAX_VALUE);
		}

	/**
	 * Create a new shared, replicated <I>bound</I> variable.
	 *
	 * @param  initialBound    Initial bound for branch-and-bound search.
	 * @param  comm            Communicator used by the <I>bound</I> variable.
	 * @param  tag             Message tag used by the <I>bound</I> variable.
	 *                         This message tag must not be used in any other
	 *                         messages in the program.
	 */
	public static ReplicatedInteger createBoundVariable
		(int initialBound,
		 Comm comm,
		 int tag)
		{
		return new ReplicatedInteger
			(IntegerOp.MINIMUM, initialBound, tag, comm);
		}

// Hidden data members.

	// List of DNA sequences with which to construct trees.
	private DnaSequenceList seqList;

	// Shared, replicated bound variable.
	private ReplicatedInteger 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, replicated <I>bound</I> variable.
	 * @param  results  Object in which to store the results.
	 */
	public MaximumParsimonyBnbHyb
		(DnaSequenceList seqList,
		 ReplicatedInteger 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.
	 * <P>
	 * The search will commence at level <I>L</I> of the search graph, 0 &le;
	 * <I>L</I> &le; <I>N</I>&minus;1, where <I>N</I> is the number of sequences
	 * in the DNA sequence list. Of the (2<I>L</I>&minus;1)!! vertices at level
	 * <I>L</I>, the search will start at the <I>V</I><SUB>1</SUB>-th such
	 * vertex and end at the <I>V</I><SUB>2</SUB>-th such vertex, 0 &le;
	 * <I>V</I><SUB>1</SUB> &le; <I>V</I><SUB>2</SUB> &le;
	 * (2<I>L</I>&minus;1)!!&nbsp;&minus;&nbsp;1. All search graph vertices at
	 * and below vertices <I>V</I><SUB>1</SUB> through <I>V</I><SUB>2</SUB>
	 * inclusive will be searched.
	 * <P>
	 * The results are stored in the {@linkplain MaximumParsimonyResults} object
	 * specified to the constructor. The <TT>findTrees()</TT> method will only
	 * find trees whose parsimony scores are less than or equal to the value of
	 * the <TT>bound</TT> variable specified to the constructor or the best
	 * bound found thereafter, whichever is smaller.
	 *
	 * @param  startLevel
	 *     <I>L</I>, the level of the search graph at which to commence the
	 *     search.
	 * @param  vertex1
	 *     <I>V</I><SUB>1</SUB>, the search graph vertex at level <I>L</I> at
	 *     which to start the search.
	 * @param  vertex2
	 *     <I>V</I><SUB>2</SUB>, the search graph vertex at level <I>L</I> at
	 *     which to end the search.
	 *
	 * @exception  IOException
	 *     Thrown if an I/O error occurred.
	 */
	public void findTrees
		(int startLevel,
		 int vertex1,
		 int vertex2)
		throws IOException
		{
		// Initialize tree signature as specified by <startLevel> and <vertex1>.
		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);

				// 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;
					}
				}
			}
		}

	}