<?php
/*
* bindist.php
* http://tasuke.dna.affrc.go.jp/
*
* Author:Norio Tabei 
*
* Copyright (c) 2024 National Institute of Agrobiological Sciences   
*
*/
/*
set_time_limit(10);
ignore_user_abort(false); 
*/
require '../conf/config.php';
require './common.php';
require './commonDB.php';

$prm_chr=$argv[1];
$prm_s=$argv[2];
$prm_e=$argv[3];
$qry_list=$argv[4];
$prm_qual=$argv[5];
$prm_dp=$argv[6];
$prm_list=$argv[7];
$prm_cri=$argv[8];
$prm_gene=$argv[9];
$prm_name=$argv[10];
$prm_obj_id=$argv[15];
$prm_phy_grp=$argv[12];
$servername=$argv[13];
$headurl=$argv[14];
$sid=$argv[11];
$variant_types=$argv[16];
$dir = __DIR__;
$list = explode(",",$prm_list);
_pregMatch($servername);
_pregMatch($headurl);
if (!isset($_SERVER["HTTP_HOST"]) && !empty($servername)) {
	$_SERVER["HTTP_HOST"] = $servername;
}
$internal_infile = "";
$internal_outfile = "";

$dn2table = array();

if (!isset($_SERVER["HTTP_HOST"])) {
    parse_str($prm_chr, $_POST);
    foreach($_POST as $key => $val){
        $_POST[$key] = rtrim($val);
    }
    $internal_infile = $_POST['tmpdir']."/".$_POST['inseq'];
    $internal_outfile = $_POST['tmpdir']."/".$_POST['outseq'];
    chdir($_POST['path']);
    set_time_limit(0);
}

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

if(isset($_SERVER["HTTP_HOST"])){
    _pregMatch($prm_obj_id);
    _pregMatch($prm_name);
    if(!_validate_chrs($conn,$prm_chr)){_stopError("||1Invalid value");}
    if(!_validate_num($prm_s)){_stopError("||2Invalid value");}
    if(!_validate_num($prm_e)){_stopError("||3Invalid value");}
    $se_len = $prm_e - $prm_s + 1;
    if($prm_cri == 1 && $prm_gene == 1 && !empty($prm_obj_id)) {
      $reg_info = getCDSLength($conn, $prm_obj_id);
      $se_len = $reg_info['cds']['len'];
    }
    $maxLength = 20000;
    if ($se_len > $maxLength) {_stopError("||4Invalid value:Length>". $maxLength);}
    if(!_validate_num($prm_qual)){_stopError("||5Invalid value");}
    if(!_validate_num($prm_dp)){_stopError("||6Invalid value");}
    if(!_validate_each_acc($list)){_stopError("||7Invalid value");}
    _pregMatchQueryList($qry_list);
    if($prm_cri != 1 and $prm_cri != 0){_stopError("||8Invalid value");}
    if($prm_gene != 1 and $prm_gene != 0){_stopError("||9Invalid value");}
    if(_chk_cds($conn, $prm_obj_id) < 1 and $prm_gene == 1){ _stopError("||This object does not have CDS region.");}
    $prm_phy_grp = _validate_group($prm_phy_grp);
    _pregMatch($sid);
}

if(!isset($_SERVER["HTTP_HOST"])){
    $seq = _convCompDelSeq($internal_infile,$prm_cri,$prm_gene,$conn,$prm_obj_id,$prm_s);
    file_put_contents($internal_outfile, $seq."\n");
    exit;
}

if($prm_cri == 1){
    $params = _get_pos_from_id($conn, $prm_obj_id);
    if(count($params) != 3){
        _stopError("||The ID was not found.");
    }
    $prm_chr = $params[0];
    $prm_s = $params[1];
    $prm_e = $params[2];
}
if(isset($qry_list)) {
	$prm_list = "";
	$qryarr = explode(",",$qry_list);
	$count = count($qryarr);
	$query = "select accession_db from sys_accession where accession_id = ?";
        $stmt = $conn->stmt_init();
        $stmt->prepare($query);
	for ($i=0; $i<$count; $i++) {
	    $stmt->bind_param("s", $qryarr[$i]);
            $stmt->execute();
            $stmt->store_result();
            $stmt->bind_result($row0);
            while ($stmt->fetch()) {
		$prm_list = $prm_list.$row0;
		if($i<$count-1) {
			$prm_list = $prm_list.",";
		}
	    }
	}
	$stmt->free_result();
	$stmt->close();
}

$wkdir = tempnam($tempdir, 'tasuke_phylotree_bindist');
unlink($wkdir);
if(mkdir($wkdir, 0700) == false) _stopError("||Can not create temporary directory.");
register_shutdown_function(function() use($wkdir) {
    // Make sure to remove working directory even if the process terminates in error.
    cleanDir($wkdir);
});

$send_obj_id = "";
if(!empty($prm_obj_id) && $prm_gene) $send_obj_id = $prm_obj_id;
$binTable = _create_infile_bindist($conn, $wkdir, $prm_chr,$prm_s,$prm_e,$prm_qual,$prm_dp,$variant_types,$prm_list,$send_obj_id);
$dmatrixPath = execMakeBinaryDistanceMatrix($wkdir, $makeBinaryDistanceMatrixPath, $binTable);
$infile = convMatrixForPhylip($wkdir, $dmatrixPath);

if(chdir($wkdir) == false){_stopError("||Can not move to temporary directory.");}

// Make Newick by NJ
$cmd = sprintf("echo -e \"L\nY\n\" | %s/exe/neighbor", $phylipDir);
exec($cmd, $resArray, $return);
//error_log("CMD:". $cmd. " => ". $return);
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);
//error_log("CMD:". $cmd. " => ". $return);
if($return !=0) {_stopError("||Can not execute PHYLIP:retree.");}

$nwk = _get_nwk($wkdir);

$order =_make_order_by_newick($conn ,$prm_name, $nwk);

$res = _make_results($order, $nwk, $prm_phy_grp, $conn, $prm_name);
echo $res;

//////////////// functions ////////////////////
function _make_results($order, $nwk, $prm_phy_grp, $conn,$name){
    $trimed_order = array();
    $relation = explode(",", $order);
    $nwk = preg_replace("/\s/", "", $nwk);
    for($i=0; $i<count($relation); $i++){
        $acc_and_grp = "";
        $accession = explode(":",$relation[$i]);
	$acc_id = $accession[0];
	$acc_name = $accession[2];
	if ($name == "name") {
        	$acc_and_grp = $acc_name."_"._get_grp($conn, $acc_name, $prm_phy_grp,$name);
	} else if ($name == "id") {
        	$acc_and_grp = $acc_id."_"._get_grp($conn, $acc_id, $prm_phy_grp,$name);
	}
        $patt = "/(\,|\)|\()$accession[1]:/";
        $patt2 = "/>$accession[1]\n/";
        array_push($trimed_order, $acc_name.":".$acc_id);
    }
    return implode(",",$trimed_order)."|".$nwk."|";
}

/*
function _clean_dir($wkdir){
    $php_tmp_dir = dirname($wkdir);
    $files = Array("infile_dnadist", "outfile_dnadist", "infile", "outfile", "outtree", "header", "phylip_wo_header", "cmd_dnadist", "cmd_tree", "cmd_retree");
    for($i=0; $i<count($files); $i++){
        unlink($files[$i]);
    }
    chdir($php_tmp_dir);
    rmdir($wkdir);
}
*/

function _get_nwk($wkdir){
    return file_get_contents($wkdir."/outtree");
}

function _make_order_by_newick($conn, $name, $newick) {
	$vals = explode(",", $newick);
	$valsLen = count($vals);
	$tailPat = '/\:[^\:]+$/';

	$order = Array();
	for($i=0; $i<$valsLen; $i++) {
		$cnvVal = $vals[$i];
		while(preg_match($tailPat, $cnvVal)) {
			$cnvVal = preg_replace($tailPat, '', $cnvVal);
		}
		$cnvVal = preg_replace('/^[\s\(]+/', '', $cnvVal);
		$cnvVal = preg_replace('/[\s\)]+$/', '', $cnvVal);
		if(preg_match('/d\d+/', $cnvVal, $matches)) {
			array_unshift($order, $matches[0]);
//error_log("Test:". $vals[$i]. " => ". $matches[0]);
		}
/*
		if(!preg_match('/^d\d+$/', $cnvVal)) continue;
		array_unshift($order, $cnvVal);
*/
	}
//error_log("Test:". $newick);

	$accessions = Array();
	for($i=0; $i<count($order); $i++){
		$str = _get_accessionid($conn, $order[$i], $name);
		array_push($accessions, $str);
	}
	return implode(",", $accessions);
}
/*
function getDb2IdMap($conn) {
	$query = "select accession_id,accession_db from sys_accession";
	$stmt = $conn->stmt_init();
	$stmt->prepare($query);
	$stmt->execute();
	$stmt->store_result();
	$stmt->bind_result($rowId, $rowDb);
	$db2id = array();
	while ($stmt->fetch()) {
		$db2id[$rowDb] = $rowId;
	}
	$stmt->free_result();
	$stmt->close();
	return $db2id;
}
*/

function convMatrixForPhylip($phylipDir, $inDmatrix) {
	$outDmatrix = $phylipDir. "/infile";

	$in = fopen($inDmatrix, 'r');
	$out = fopen($outDmatrix, 'w');

	$line = fgets($in);
	fwrite($out, $line);	// matrix header(sample num)
	while($line = fgets($in)) {
		$line = rtrim($line);
		if(!preg_match("/\S/", $line) || preg_match("/^[\s]+$/", $line)) continue;

		$onerowArray = preg_split("/\s+/", $line);
		$onerow = "";
		if(count($onerowArray) < 2) {
			$onerow = sprintf("%-9s \n", $onerowArray[0]);
		} else {
			$smplName = array_shift($onerowArray);
			$onerow = sprintf("%-9s %s\n", $smplName, implode(" ", $onerowArray));
		}
		fwrite($out, $onerow);
	}

	fclose($out);
	fclose($in);

	return $outDmatrix;
}
function execMakeBinaryDistanceMatrix($wkdir, $binPath, $binTable) {
	$dmatrixPath = $wkdir. "/dmatrix";
	$error_log = $wkdir. "/stderr";

	$cmd = $binPath. " -i ". $binTable. " >". $dmatrixPath. " 2>". $error_log;
	exec($cmd, $resArray, $return);
	if($return !=0) _stopError("||Can not execute makeBinaryDistanceMatrix.");

	return $dmatrixPath;
}
function getSqlParams($preParams, $isPlus, $query, $prm_chr, $prm_qual, $prm_dp, $variant_types_str) {
	if($prm_qual == 0 || !empty($prm_qual)) {
		if($isPlus) {
			$query .= " and QUAL>=?";
		} else {
			$query .= " and GQ>=?";
		}
		$preParams[0] .= "i";
		$preParams[] = $prm_qual;
	}
	if($prm_dp == 0 || !empty($prm_dp)) {
		$query .= " and dp>=?";
		$preParams[0] .= "i";
		$preParams[] = $prm_dp;
	}
	if(!empty($variant_types_str)) {
		$typesSql = "";
		$hitSid = 0;
		if(preg_match("/s/i", $variant_types_str)) {
			$typesSql .= "'s',";
			$hitSid++;
		}
		if(preg_match("/i/i", $variant_types_str)) {
			$typesSql .= "'i',";
			$hitSid++;
		}
		if(preg_match("/d/i", $variant_types_str)) {
			$typesSql .= "'d',";
			$hitSid++;
		}
		$typesSql =  preg_replace("/,$/", "", $typesSql);

		// if($variant_types_str == 'sid'), Since all types are targeted, there is no need to filter.
		if(!empty($typesSql) && $hitSid>=1 && $hitSid<=2) $query .= " and variant in (". $typesSql. ")";
	}
//error_log("Test:". $query. " <= ". implode(",", $preParams));
	return array($query, $preParams);	
}
function _create_infile_bindist($conn, $wkdir, $prm_chr,$prm_s,$prm_e,$prm_qual,$prm_dp,$variant_types,$prm_list, $prm_obj_id){
    $list=explode(",", $prm_list);
    $enabledEffImpacts = array("0" => True, "1" => True, "2" => True, "3" => True, "-1" => True);
    $uniqPOS = _allLIST($conn, $list, $prm_qual, $prm_dp, $variant_types, $prm_chr, $prm_s, $prm_e, $enabledEffImpacts, $prm_obj_id);

    $isPlus = modeIsPlus($conn);

    $binTable = $wkdir. "/bintable";
    $fp = fopen($binTable, "w");

//$db2id = getDb2IdMap($conn);
    $numUniPOS = count($uniqPOS);
    $list_cnt = count($list);

    // Output Reference row
    fwrite($fp, "Reference". str_repeat(",0", $numUniPOS). "\n");

    for ( $i=0; $i<$list_cnt; $i++ ) {
	fwrite($fp, $list[$i]);

	$stmt = $conn->stmt_init();
	//TODO There is a way to speed it up (same as DB-scanning of tasuke_tree_dmatrix.pl, same as multi-shread of get_fasta_multi.php)
	$query = "";
	if($isPlus) {
		$query = "SELECT GT FROM ". $list[$i]. "_variant WHERE ". getGtNotZeroSql(). " and chr=? and pos=?";
	} else {
		//TODO Is sys_variant join necessary?
		$query = "SELECT GT FROM sys_variant s,". $list[$i]. "_variant v WHERE ". getGtNotZeroSql().
			" and s.chr=v.chr and s.pos_vcf=v.pos_vcf and v.chr=? and v.pos_vcf=?";
	}
	list($query, $preParams) = getSqlParams(array("si", $prm_chr, NULL), $isPlus, $query, $prm_chr, $prm_qual, $prm_dp, $variant_types);
	$params = array();
	$numSqlParams = count($preParams);
	for($j=0; $j<$numSqlParams; $j++) { $params[$j] = &$preParams[$j]; }
	$stat = $stmt->prepare($query);
	for ( $j=0; $j<$numUniPOS; $j++ ) {
		$params[2] = &$uniqPOS[$j];
		call_user_func_array(array($stmt, 'bind_param'), $params);
//		$stmt->bind_param("si", $prm_chr, $uniqPOS[$j]);
		$stmt->execute();
		$stmt->store_result();
		$stmt->bind_result($rowGt);
		$bin;
		if($stmt->fetch()) {
			if($rowGt === "0/0" || empty($rowGt)) {
				$bin = "0";
			} else if($rowGt === "./.") {
				$bin = "NA";
			} else {
				$chkDPStr = _CHKDP($list[$i], $prm_chr, $uniqPOS[$j], $conn);
				if ( empty($chkDPStr) || !empty($rowGt)) {
					// DEPTH != NA
					if($rowGt === "1/1" || $rowGt === "1/2") {
						$bin = "1";
					} else if($rowGt === "0/1") {
						$bin = "2";
					} else {
						$bin = "0";
					}	
				} else {
					// DEPTH = NA
					$bin = "NA";
				}
			}
		} else {
			$bin = 0;
		}
//error_log(($j+1). "/". $numUniPOS. " ". $db2id[$list[$i]]. " ". $prm_chr. " ". $uniqPOS[$j]. " ". $rowGt. " => ". $bin);
		fwrite($fp, ",". $bin);
		$stmt->free_result();
	}
	fwrite($fp, "\n");
	$stmt->close();
    }
    fclose($fp);

    return $binTable;
}
function _allLIST($conn, $list, $prm_qual, $prm_dp, $variant_types, $prm_chr, $prm_s, $prm_e, $enabledEffImpacts, $prm_obj_id){
	$isPlus = modeIsPlus($conn);

	$stmt = $conn->stmt_init();
	$numList = count($list);
	$uniPos = array();
	for($i=0; $i<$numList; $i++){
		if($isPlus) {
			$query="SELECT pos,EFF,QUAL,DP FROM ". $list[$i]. "_variant WHERE ". getGtNotZeroSql(). " and chr=? and pos>=? and pos<=?";
		} else {
			$query="SELECT v.pos_vcf,EFF,GQ,DP FROM sys_variant s,". $list[$i]. "_variant v WHERE ".
				getGtNotZeroSql(). " and s.chr=v.chr and s.pos_vcf=v.pos_vcf and v.chr=? and v.pos>=? and v.pos<=?";
		}
//$query .= " and pos not in(18621,18622,18693,18694,18699,18700)";
		//TODO VarianType filter
		list($query, $preParams) = getSqlParams(array("sii", $prm_chr, NULL, NULL), $isPlus, $query, $prm_chr, $prm_qual, $prm_dp, $variant_types);
		$params = array();
		$numSqlParams = count($preParams);
		for($j=0; $j<$numSqlParams; $j++) { $params[$j] = &$preParams[$j]; }
		$stmt->prepare($query);
		$params[2] = &$prm_s;
		$params[3] = &$prm_e;
		call_user_func_array(array($stmt, 'bind_param'), $params);
//		$stmt->bind_param("sii", $prm_chr, $prm_s, $prm_e);
		$stmt->execute();
		$stmt->store_result();
		$stmt->bind_result($rowPos, $rowEff, $rowQual, $rowDp);
//$partPos = array();
		while ($stmt->fetch()) {
			$uniPos[strval($rowPos)] = True;
//$partPos[strval($rowPos)] = True;
		}
//error_log("Test:". $list[$i]. ",". count($partPos). "/". count($uniPos));
		$stmt->free_result();
	}
	$stmt->close();

	$preResArr = array_keys($uniPos);
	sort($preResArr);
	$resArr = array();

	// CDS region filter
	if($prm_obj_id) {
		$cdsRegs = _get_cds_region($conn, $prm_obj_id);
		$numCdsRegs = count($cdsRegs);
		$numPreResArr = count($preResArr);
		for($i=0; $i<$numPreResArr; $i++) {
			$inCds = False;
			for($j=0; $j<$numCdsRegs; $j++) {
				if($preResArr[$i] >= $cdsRegs[$j][0] && $preResArr[$i] <= $cdsRegs[$j][1]) {
					$inCds = True;
					break;
				}
			}
			if($inCds) $resArr[] = $preResArr[$i];
//if(!$inCds) error_log("SkipPos:". $preResArr[$i]);
		}
	} else {
		$resArr = $preResArr;	// It is not a copy in memory(PHP's Copy on Write=pass by reference unless it is changed).
	}
	return $resArr;
}

function _chk_cds($conn, $prm_obj_id){
    $cds = 0;
    $query = "select seq from sys_track where id=? and (type='CDS' or type='cds') limit 1";
    $stmt = $conn->stmt_init();
    $stmt->prepare($query);
    $stmt->bind_param("s", $prm_obj_id); 
    $stmt->execute();
    $stmt->store_result();
    $stmt->bind_result($row0);
    while ($stmt->fetch()) {
        $cds++;
    }
    $stmt->free_result();
    $stmt->close();
    return $cds;
}

function _get_cds_region($conn, $prm_obj_id){
    $cdsRegs = array();
    $query = "select pos,end from sys_track where id=? and (type='CDS' or type='cds') order by pos";
    $stmt = $conn->stmt_init();
    $stmt->prepare($query);
    $stmt->bind_param("s", $prm_obj_id); 
    $stmt->execute();
    $stmt->store_result();
    $stmt->bind_result($row0,$row1);
    while ($stmt->fetch()) {
        $cdsRegs[] = array($row0, $row1);
    }
    $stmt->free_result();
    $stmt->close();
    return $cdsRegs;
}

function _get_pos_from_id($conn, $prm_obj_id){
    $pos = array();
    $query = "select chr,pos,end from sys_track where id=? and parent=1 limit 1";
    $stmt = $conn->stmt_init();
    $stmt->prepare($query);
    $stmt->bind_param("s", $prm_obj_id); 
    $stmt->execute();
    $stmt->store_result();
    $stmt->bind_result($row0,$row1,$row2);
    while ($stmt->fetch()) {
        $pos[0] = $row0;
        $pos[1] = $row1;
        $pos[2] = $row2;
    }
    $stmt->free_result();
    $stmt->close();
    return $pos;
}

function _validate_chrs($conn, $prm_Chr){
    $flag = false;
    $query = "select seq from sys_chromosome where id=? limit 1";
    $stmt = $conn->stmt_init();
    $stmt->prepare($query);
    $stmt->bind_param("s", $prm_Chr); 
    $stmt->execute();
    $stmt->store_result();
    $stmt->bind_result($row0);
    while ($stmt->fetch()) {
        $flag = true;
    }
    $stmt->free_result();
    $stmt->close();
    if($flag){
       return true;
    }else{
       return false;
    }
}

function _validate_each_acc($list){
    $num = count($list);
    for($a=0;$a<$num;$a++){
        if (preg_match("/^d\d+$/i", $list[$a])) {
            //pass
        }else{
           return false;
        }
    }
    return true;
}

function _validate_each_qry($list){
    $num = count($list);
    for($a=0;$a<$num;$a++){
        if (preg_match("/\;|\||\&|\`|\(|\)|\<|\>|\*|\?|\{|\}|\[|\]|\!/", $list[$a])) {
            //pass
        }else{
           return false;
        }
    }
    return true;
}

function _validate_name($str){
    if($str == "id"){
        $str="accession_id";
    }else if($str == "name"){
        $str="accession_name";
    }else{
        return false;
    }
    return true;
}

function _validate_group($str){
    $keys = array(
        "variety" => "1",
        "sub_variety" => "1",
        "origin" => "1",
        "origin_2" => "1",
        "type" => "1",
        "other_1" => "1",
        "other_2" => "1",
    );
    if($str == "subvariety" ){$str = "sub_variety";}
    if(array_key_exists($str, $keys)){
        return $str;
    }else{
        return false;
    }
}

function _get_accessionid($conn, $internal_id,$name){
    $id = "";
    $query = "select accession_id,accession_name from sys_accession where accession_db=? limit 1";
    $stmt = $conn->stmt_init();
    $stmt->prepare($query);
    $stmt->bind_param("s", $internal_id); 
    $stmt->execute();
    $stmt->store_result();
    $stmt->bind_result($row0,$row1);
    while ($stmt->fetch()) {
        $id = $row0;
	$name = $row1;
    }
    $stmt->free_result();
    $stmt->close();
    //return substr($id, 0, 10).":".$internal_id;
    return $id.":".$internal_id.":".$name;
}


function _get_grp($conn, $id, $grp,$name){
    $ret = "";
    if ( $name == "id") {
    	$query = "select ".$grp." from sys_accession where accession_id=? limit 1";
    } else if($name == "name") {
    	$query = "select ".$grp." from sys_accession where accession_name=? limit 1";
    }
    $stmt = $conn->stmt_init();
    $stmt->prepare($query);
    $stmt->bind_param("s", $id); 
    $stmt->execute();
    $stmt->store_result();
    $stmt->bind_result($row0);
    while ($stmt->fetch()) {
        $ret = $row0;
    }
    $stmt->free_result();
    $stmt->close();
    return $ret;
}
function _CHKDP($a, $b, $c, $conn){
    global $db, $dn2table;

    $table = "";
    if(array_key_exists($a, $dn2table)) {
        $table = $dn2table[$a];
    } else {
        $table = getSerialTableByColumn($conn, $db, $a."_accession", $b."_b");
        $dn2table[$a] = $table;
    }
    $query="SELECT pos FROM ". $table. " WHERE `".$b."_b` is null and pos= ? ";
//    $query="SELECT pos FROM ".$a."_accession WHERE `".$b."_b` is null and pos= ? ";
    $stmt = $conn->stmt_init();
    $stat = $stmt->prepare($query);
    if(!$stat) { return NULL; }
    $stmt->bind_param("i", $c);
    $stmt->execute();
    $stmt->store_result();
    $stmt->bind_result($row0);
    $stmt->fetch();
    $rtn = $row0;
    $stmt->free_result();
    $stmt->close();
    return $rtn;
}

function returnError($str){
    echo $str;
}

function _stopError($str){
    echo "Error: ". $str;
    exit;
}

?>
