<?php
/*
* getNewick.php
* http://tasuke.dna.affrc.go.jp/
*
* Author: Norio Tabei
*
* Copyright (c) 2020 National Institute of Agrobiological Sciences
*
*/
//ini_set("memory_limit", "512M");

require '../conf/config.php';
require './common.php';
require_once "./PHPPhylogeneticTrees/autoload.php";
require 'globalVars.php';

use PHPPhylogeneticTrees\NJTree;

if(!file_exists($distanceMatrixPath)) {
	echo "{}";
	exit;
}

$options = getopt('d');
if($options === false){
	$options = array();
}

//////////////////////////////////  Connecting with DB
if( !$conn = new mysqli($host, $user, $pswd, $db)) {
    error_log("Connection failed.\n");
    exit;
}

$accDbs;
(isset($_POST['prm_Lists']) && !empty($_POST['prm_Lists'])) ? $accDbs = explode(",", $_POST['prm_Lists']) : $accDbs = getAllAccDbs($conn);
//(isset($_POST['prm_Lists']) && !empty($_POST['prm_Lists'])) ? $prm_Lists = $_POST['prm_Lists'] : $prm_Lists = implode(",", getAllAccDbs($conn));
//(isset($_POST['prm_Lists']) && !empty($_POST['prm_Lists'])) ? $prm_Lists=$_POST['prm_Lists'] : returnError("Error");

$idConvMapArr = getId2DbMap($conn);
$id2DbMap = $idConvMapArr[0];
$db2IdMap = $idConvMapArr[1];

$orders = getOrders();
$numOrders = count($orders);
$orderAccIds = array();
foreach($accDbs as $accdb) {
	if(!array_key_exists($accdb, $db2IdMap)) continue;
	$accid = $db2IdMap[$accdb];
	if($numOrders === 0 || array_key_exists($accid, $orders)) $orderAccIds[] = $accid;
}

$taxa_dmatrix = readDistanceMatrixFile($distanceMatrixPath);
$taxa = $taxa_dmatrix[0];
$d_mat = $taxa_dmatrix[1];

// Holds a callback to pass to NJTree.php
abstract class AbstractCallback {
	abstract public function call(); 
}

if(count($orderAccIds) > 0) {
	$res = rebuildMatrix($taxa, $d_mat, $orderAccIds);
	$taxa = $res[0];
	$d_mat = $res[1];
}

//var_dump($taxa);
//var_dump($d_mat);

$taxaDb = array();
foreach($taxa as $id) {
	$taxaDb[] = $id2DbMap[$id];
}
if(!isset($systreeProgressChar)) $systreeProgressChar = DEF_SYSTREE_PROGRESS_CHAR;
if(!isset($systreeProgressCharLen)) $systreeProgressCharLen = DEF_SYSTREE_PROGRESS_CHAR_LEN;

$newick = "";
$needrooting = false;
if($systreeUsePhylip !== 0 && !empty($phylipDir) && file_exists($phylipDir. "/exe/neighbor")) {
	// use Phylip neighbor.
	$tmppath = tempnam($tempdir, 'tasuke-');
	unlink($tmppath);

	$wkdir = $tmppath;
//error_log("wkdir:". $wkdir. "  matrix:". $distanceMatrixPath);
	if(mkdir($wkdir) == false){_stopError("||Can not create temporary directory.");}
	if(chdir($wkdir) == false){_stopError("||Can not move to temporary directory.");}
	register_shutdown_function(function() use($wkdir) {
//		_clean_dir($wkdir);
		cleanDir($wkdir);
	});
	$infile = $wkdir. "/infile";
	createPhylipInfile($taxaDb, $d_mat, $infile);

	$progUnitStr = str_repeat($systreeProgressChar, $systreeProgressCharLen);

	$cmd = sprintf("echo -e \"%s\nL\nY\n\" | %s/exe/neighbor", $infile, $phylipDir);
	$handle = popen($cmd, "r");
	while (!feof($handle)) {
		$buffer = fgets($handle, 4096);
		if(preg_match('/^Cycle/i', $buffer)) {
			if(array_key_exists("d", $options)) {
				fputs(STDERR, $buffer);
			} else {
				echo $progUnitStr;
			}
			if(ob_get_level() > 0) ob_flush();
			flush();
		}
	}
	$return = pclose($handle);
	if($return !=0) {_stopError("||Can not execute PHYLIP:neighbor.");}

	// Mid-point rooting
	$cmd = sprintf("echo -e \"Y\n%s/outtree\nM\nW\nR\nR\nQ\n\" | %s/exe/retree > /dev/null", $wkdir, $phylipDir);
	exec($cmd, $resArray, $return);
	if($return !=0) {_stopError("||Can not execute PHYLIP:neighbor.");}

	$nwk = file_get_contents($wkdir. "/outtree");
	$nwk = preg_replace("/\s/", "", $nwk);
//error_log($nwk);

	$newick = $nwk;
	$needrooting = false;
/*
	if(array_key_exists("d", $options)) {
		echo $nwk;
	} else {
		echo json_encode(array("newick" => $nwk, "needrooting" => false));
	}
*/
} else {
	// use PHPPhylogeneticTrees module.
	$numDmat = count($d_mat);
	for($i=0; $i<$numDmat; $i++) {
		$d_mat[$i] = $d_mat[$i]->toArray();	// Convert SplFixedArray[] back to a 2D-array() to set to NJTree. It works without this, but just in case.
	}

	// If you want NJTree to do additional processing, pass a class that contains callbacks instead of writing directly.
	// Anonymous functions are easier, but don't use them as they are supported since PHP 5.3
	$callbacks = array();
	class CallbackCycle extends AbstractCallback {
		private $strUnit = "";
		private $opt = array();
		public function __construct($str, $times, $opt) {
			$this->strUnit = str_repeat($str, $times);
			$this->opt = $opt;
		}
		public function call() {
			if(!array_key_exists("d", $this->opt)) {
				echo $this->strUnit;
			}
			if(ob_get_level() > 0) ob_flush();
			flush();
		}
	}
	$callbacks['cycle'] = new CallbackCycle($systreeProgressChar, $systreeProgressCharLen, $options);
	$nj = new NJTree($d_mat, $callbacks);
//error_log("Test2:". json_encode($nj->getLabels()));
	$newick = $nj->getNewick($taxaDb);

	// Mid-point rooting is done on the browser side.
	$needrooting = true;
/*
	if(array_key_exists("d", $options)) {
		echo $newick;
	} else {
		echo json_encode(array("newick" => $newick, "needrooting" => true));
	}
*/
}
//error_log("Max:". memory_get_peak_usage());

if(array_key_exists("d", $options)) {
	// Since the phylip format cannot handle sample names with more than 10 characters,
	// create a Newick and then replace AccDB with AccID.
	$dbOrder = $taxaDb;	// Deep copy
	array_multisort( array_map( "strlen", $dbOrder), SORT_DESC, $dbOrder);	// Destructive sort
	foreach($dbOrder as $db) {
		$id = $db2IdMap[$db];
//		echo "$db => $id\n";
		$newick = preg_replace("/". $db. "/", $id, $newick);
	}
	echo $newick;
} else {
	echo json_encode(array("newick" => $newick, "needrooting" => $needrooting));
}
function getOrders() {
	$orders = array();
	if (file_exists( "../conf/order.conf")) {
		$line = file_get_contents("../conf/order.conf");
		$accids = explode("\n",$line);
		foreach ($accids as $accid) {
			$accid = trim($accid);
			if(empty($accid)) {continue;}
			$orders[$accid] = 1;
		}
	}
	return $orders;
}
function getAllAccDbs($conn) {
	$dbAccs = array();
	$query = "SELECT accession_id,accession_db FROM sys_accession";
	$stmt = $conn->stmt_init();
	$stmt->prepare($query);
	$stmt->execute();
	$stmt->store_result();
	$stmt->bind_result($accid, $accdb);
	while ($stmt->fetch()) {
		$dbAccs[] = $accdb;
/*
		if($numOrders === 0) {
			$dbAccs[] = $accdb;
		} else if($numOrders > 0 && array_key_exists($accid, $orders)) {
			$dbAccs[] = $accdb;
		}
*/
	}
	return $dbAccs;
}
function getId2DbMap($conn) {
	$id2db = array();
	$db2id = array();
	$stmt = $conn->stmt_init();
	$query = "SELECT accession_id, accession_db FROM sys_accession";
	$stmt->prepare($query);
	$stmt->execute();
	$stmt->store_result();
	$stmt->bind_result($id, $db);
	while ($stmt->fetch()) {
		$id2db[$id] = $db;
		$db2id[$db] = $id;
	}
	return array($id2db, $db2id);
}
function readDistanceMatrixFile($file) {
//error_log("start:". memory_get_usage());
	$taxa = array();
	$dmatrix = array();
	$IN = fopen($file, "r");
	$tgtIndex = 0;
	if($IN){
		while($line = fgets($IN)) {
			$line = rtrim($line);
			if(!preg_match("/\S/", $line) || preg_match("/^[\s\d]+$/", $line)) continue;

			$onerowArrayPre = preg_split("/\s+/", $line);	// preg_split("/\s+/") is heavier than explode(" "), but is required for strict Phylip formatting.
			$onerowArray = new SplFixedArray($tgtIndex);	// Use SplFixedArray because PHP's 2D-array() uses too much memory.
			$acc = $onerowArrayPre[0];
			$taxa[] = $acc;
			for($i=1; $i<=$tgtIndex; $i++) {
				$onerowArray[$i-1] = (float) $onerowArrayPre[$i];
			}
			$dmatrix[] = $onerowArray;
			$onerowArrayPre = null;
			unset($onerowArrayPre);
			$tgtIndex++;

/*
			$onerowArrayPre = preg_split("/\s+/", $line);
			$onerowArray = [];
			$acc = array_shift($onerowArrayPre);
			$taxa[] = $acc;
			$arrayLen = count($onerowArrayPre);
			for($i=0; $i<$tgtIndex && $i<$arrayLen; $i++) {
				$onerowArray[$i] = (float) $onerowArrayPre[$i];
			}
			$dmatrix[] = $onerowArray;
			$tgtIndex++;
*/
		}
	}
	fclose($IN);
//error_log("end:". memory_get_usage());
	return array($taxa, $dmatrix);
}
function rebuildMatrix($taxa, $d_mat, $orderAccIds) {
	if(count($orderAccIds) === 0) return array($taxa, $d_mat);

	$taxaLen = count($taxa);
	$taxaId2Index = array();
	for($i=0; $i<$taxaLen; $i++) {
		$taxaId2Index[$taxa[$i]] = $i;
	}

	// Remove Acc that does not exist in distance matrix
	$orderAccIdsLen = count($orderAccIds);
	$destTaxa = array();
	for($i=0; $i<$orderAccIdsLen; $i++) {
		$accid= $orderAccIds[$i];
		if(array_key_exists($accid, $taxaId2Index)) $destTaxa[] = $accid;
	}
/*
	$preDestTaxa = explode(",", $prm_Lists);
	$preDestTaxaLen = count($preDestTaxa);
	$destTaxa = array();
	for($i=0; $i<$preDestTaxaLen; $i++) {
		$accid = $db2IdMap[$preDestTaxa[$i]];
		if(array_key_exists($accid, $taxaId2Index)) $destTaxa[] = $accid;
	}
*/

	$destTaxaLen = count($destTaxa);
	$destDmat = array();
	for($i=0; $i<$destTaxaLen; $i++) {
		$taxaRowIndex = $taxaId2Index[$destTaxa[$i]];
		$elem = new SplFixedArray($i);
		for($j=0; $j<$i; $j++) {
			$taxaColIndex = $taxaId2Index[$destTaxa[$j]];
			$d;
			if($taxaRowIndex < $taxaColIndex) {
				$d = $d_mat[$taxaColIndex][$taxaRowIndex];
			} else {
				$d = $d_mat[$taxaRowIndex][$taxaColIndex];
			} 
//echo "$i => $taxaRowIndex, $j => $taxaColIndex, = $d\n";
			$elem[$j] = $d;
		}
		$destDmat[$i] = $elem;
	}	
	return array($destTaxa, $destDmat);
}
function createPhylipInfile($taxa, $d_mat, $infileName) {
	$fp = fopen($infileName, 'w');
	$taxaCnt = count($taxa);
	fwrite($fp, $taxaCnt."\n"); 
	for($i=0; $i<$taxaCnt; $i++) {
		$onerow = sprintf("%-9s %s\n", $taxa[$i], implode(" ", $d_mat[$i]->toArray()));
//error_log($onerow);
		fwrite($fp, $onerow);
	}
	fclose($fp);
}
function _clean_dir($wkdir){
	$php_tmp_dir = dirname($wkdir);
	chdir($wkdir);
	$files = Array("infile", "outfile", "outtree");
	for($i=0; $i<count($files); $i++){
		unlink($files[$i]);
	}
	chdir($php_tmp_dir);
	rmdir($wkdir);
}
function returnError($str){
	echo $str;
}
function _stopError($str){
	returnError("Error: ". $str);
	exit;
}
?>
