#
# tasuke_tree_dmatrix.pl
# http://tasuke.dna.affrc.go.jp/
#
# Author: Norio Tabei
#
# Copyright (c) 2021 National Institute of Agrobiological Sciences
#
package TasukeVariantExistMap;

use strict;
use Cwd;
use Data::Dumper;
use POSIX qw(ceil);
use FindBin qw($Bin);
use File::Temp qw/ tempfile /;
use File::Spec;
use Fcntl;
use NDBM_File;
use SDBM_File;
#use GDBM_File;
use lib "$Bin";
use CommonLibTASUKE;

#use Devel::Size qw(size total_size);	# This module is not installed by default, so comment it out.

sub new {
	my $class = shift;
	my $self = {
		'chrs' => [],
		'accs' => [],
		'min_qual' => 20,
		'min_dp' => 2,
		'min_dp4' => 2,
		'check_depth_null' => 0,
		'depth_null_in_memory' => 0,
		'accdb2col2table' => {},
		'verbose' => 0,
		'tempdir' => '/tmp',
		'targetv' => {'s' => 1, 'i' => 1, 'd' => 1},
		@_
	};
	return bless $self, $class;
}
sub verbose {
	my $self = shift;
	my($flg) = @_;
	$self->{'verbose'} = $flg;
}
sub setTempDir {
	my $self = shift;
	my($tempdir) = @_;
	if($tempdir) {
		if(-d $tempdir) { 
			$self->{'tempdir'} = $tempdir;
			printf STDERR "TempDir: %s\n", $tempdir  if($self->{'verbose'});
		} else {
			printf STDERR "TempDir: %s not found, using \"%s\" instead.\n", $tempdir, $self->{'tempdir'};
		}
	}
}
sub setChrs {
	my $self = shift;
	my($chrs) = @_;
	printf STDERR "Chrs: %s\n", join(",", @{$chrs})  if($self->{'verbose'});
	$self->{'chrs'} = $chrs;
}
sub setChr {
	my $self = shift;
	my($chrId) = @_;
	$self->setChrs([$chrId]);
}
sub setAccs {
	my $self = shift;
	my($accs) = @_;
	$self->{'accs'} = $accs;
}
sub setAccsByOrderFile {
	my $self = shift;
	my($file) = @_;

	printf STDERR "OrderFile: %s\n", $file  if($self->{'verbose'});
	if( !$file || !(-e $file)) {
		print STDERR "File not fount: $file\n";
		exit(1);
	}

	my @accs;
	local *FH;
	open(FH, "<$file");
	while(<FH>) {
		chomp;
		my $onerow = $_;
		next if($onerow !~ m/\S/);

		$onerow =~ s/^\s+//; 
		my @vals = split /\s+/, $onerow;
		my $acc = $vals[0];
		$acc =~ s/\W+$//;
		push @accs, $acc if($acc);
	}
	close(FH);
	
	$self->setAccs(\@accs);
	if($self->{'verbose'}) {
		print STDERR "Target accessions:\n";
		my $accsNum = @accs;
		for(my $i=0; $i<$accsNum; $i++) {
			printf STDERR "%d\t%s\n", $i+1, $accs[$i];
		}
		print STDERR "\n";
	}
}
sub setMinQual {
	my $self = shift;
	my($val) = @_;
	$self->{'min_qual'} = $val;
}
sub setMinDp {
	my $self = shift;
	my($val) = @_;
	$self->{'min_dp'} = $val;
}
sub setMinDp4 {
	my $self = shift;
	my($val) = @_;
	$self->{'min_dp4'} = $val;
}
sub setCheckDepth {
	my $self = shift;
	my($val) = @_;
	printf STDERR "CheckDepth: %s\n", $val  if($self->{'verbose'});
	$self->{'check_depth_null'} = $val;
}
sub setDepthInMem {
	my $self = shift;
	my($val) = @_;
	printf STDERR "Put Depth=Null map in memory: %s\n", $val  if($self->{'verbose'});
	$self->{'depth_null_in_memory'} = $val;
}
sub replaceTargetVariantTypes {
	my $self = shift;
	my ($hashrefTypes) = @_;	# ex. {'s' => 1, 'i' => 0, 'd' => 0}
	$self->{'targetv'} = $hashrefTypes;

	my @types;
	push @types, "SNP" if($self->{'targetv'}->{'s'});
	push @types, "INS" if($self->{'targetv'}->{'i'});
	push @types, "DEL" if($self->{'targetv'}->{'d'});
	printf STDERR "Target variant type: %s\n", join(',', @types);
}
sub replaceTargetVariantTypesFromStr {
	my $self = shift;
	my ($vtypesStr) = @_;	# 'si', 'd', etc

	$vtypesStr = lc($vtypesStr);

	my $err = 0;

	if(!$vtypesStr) {
		return 1;
	}

	my @vtypeArr = split //, $vtypesStr;
	my $numVtypeArr = @vtypeArr;
	if($numVtypeArr > 3 || $numVtypeArr < 1) {
		$err++;
	} else {
		my $hashrefVtypes = {'s' => 0, 'i' => 0, 'd' => 0};
		my $numVtypes = 0;
		foreach my $vtype (@vtypeArr) {
#			if($vtype ne 's' && $vtype ne 'i' && $vtype ne 'd') {
			if( !exists($hashrefVtypes->{$vtype}) ) {
				$err++;
				last;
			}
			$hashrefVtypes->{$vtype} = 1;
			$numVtypes++;
		}
		if($numVtypes < 1) {
			$err++;
		} else {
			$self->replaceTargetVariantTypes($hashrefVtypes);
		}
	}
	return $err;
}
sub setTargetVariantType {
	my $self = shift;
	my ($type) = @_;
	return if($type ne 's' && $type ne 'i' && $type ne 'd');
	$self->{'targetv'}->{$type} = 1;
}
sub delTargetVariantType {
	my $self = shift;
	my ($type) = @_;
	return if($type ne 's' && $type ne 'i' && $type ne 'd');
	$self->{'targetv'}->{$type} = 0;
}
sub _getFileHandle {
    my($VCF) = @_;
    my $fh;
    if ($VCF =~ m/\.tar\.gz$/) {
        open $fh, "tar zxf $VCF -O 2> /dev/null|" or die "file open error: $!";
    } elsif($VCF =~ m/\.tar\.bz2$/) {
        open $fh, "tar jxf $VCF -O 2> /dev/null|" or die "file open error: $!";
    } elsif($VCF =~ m/\.gz$/) {
        open $fh, "gzip -dc $VCF 2> /dev/null|" or die "file open error: $!";
    } elsif($VCF =~ m/\.bz2$/) {
        open $fh, "bzip2 -dc $VCF 2> /dev/null|" or die "file open error: $!";
    } else {
        open $fh, '<', $VCF or die "file open error: $!";
    }
    return $fh;
}
sub getSamples {
	my $self = shift;
	my($vcf) = @_;

	my @smpls;
	my $in = _getFileHandle($vcf);
	while(<$in>) {
		chomp;
		my $datarow = $_;
		next if($datarow !~ m/\S/);

		if($datarow =~ m/^#CHROM/) {
			my @vals = split /\t/, $datarow;
			@smpls = @vals[9 .. $#vals];
			last;
		} elsif($datarow =~ m/^#/) {
			# Nothing to do.
		} else {
			# "#CHROM" line does not exist
			last;
		}
	}
	close($in);

	return \@smpls;
}
sub getGtBin {
	my($gt) = @_;

	if($gt !~ m/[^0\|\/]/) {
#	if($gt eq '0/0' || $gt eq '0|0') {
		return '0';
	}

	my @gtParam = split(/[\/\|]/, $gt);
	return 'NA' if($gtParam[0] eq '.');

	if(@gtParam > 2) {
		my %tmphash;
		@tmphash{@gtParam} = undef;	# Slice: ArrayVals => UniqKeys
		@gtParam = sort {$a <=> $b} keys(%tmphash);
		push(@gtParam, $gtParam[0]) if(@gtParam == 1);
	}

	if($gtParam[0] eq '0' || (@gtParam > 1 && $gtParam[1] eq '0')) {
		return '2';	# 0/1:Hetero
	} else {
		return '1';	# 1/1:Alt Homo, 1/2:Alt Hetero(without Ref)
	}
}
sub getGtBinFilterVariantType {
	my $self = shift;
	my($ref, $alts, $gt) = @_;

	my $bin = getGtBin($gt);
	if($bin eq '1' || $bin eq '2') {
		# Has variant. Check varaint type filter.
		my @gtParam = split(/[\/\|]/, $gt);
		$gtParam[0] = int($gtParam[0]);
		$gtParam[1] = int($gtParam[1]);
#		my @alts = split /,/, $allAlts;
		my $refLen = length($ref);
		my $altLen;
		if($bin == 1) {
			# Alt homo
			$altLen = length($alts->[$gtParam[0] - 1]);
		} else {
			# Alt Hetero
			if($gtParam[0] != 0 && $gtParam[1] != 0) {
				# GT:1/2
				$altLen = length($gtParam[0] < $gtParam[1] ? $alts->[$gtParam[0] - 1] : $alts->[$gtParam[1] - 1]);
			} else {
				# GT:0/1
				$altLen = length($gtParam[0] == 0 ? $alts->[$gtParam[1] - 1] : $alts->[$gtParam[0] - 1]);
			} 
		}
		if(	$refLen == $altLen && !$self->{'targetv'}->{'s'} ||
			$refLen < $altLen && !$self->{'targetv'}->{'i'} ||
			$refLen > $altLen && !$self->{'targetv'}->{'d'}
		) {
			$bin = 0;
		}
	} 
	return $bin;
}
sub getVcfGtIndex {
	my($format) = @_;
	my $gtIdx = -1;
	if($format =~ m/^GT:/i) {
		$gtIdx = 0;
	} else {
		my @formatVals = split /\:/, $format;
		my $i = 0;
		foreach my $formatVal (@formatVals) {
			if($formatVal eq 'GT') {
				$gtIdx = $i;
				last;
			}
			$i++;
		}
	}
	return $gtIdx;
}
sub getBinBitMaps {
	my %bin2bit = (
		'0' => '00',
		'1' => '01',
		'2' => '10',
		'NA' => '11'
	);
	my %bit2bin;
	foreach my $bin (keys(%bin2bit)) {
		$bit2bin{$bin2bit{$bin}} = $bin;
	}

	return (\%bin2bit, \%bit2bin);
}

###
 # Standard methods for multi-sumple VCF file.
###
sub outputFromMvcf {
	my $self = shift;
	my($vcf, $out, $accOrder) = @_;

	$accOrder = $self->getSamples($vcf) if(!$accOrder);	# Include null element.
	my @validAccs;
	my $numAccOrder = @{$accOrder};
	for(my $i=0; $i<$numAccOrder; $i++) {
		push @validAccs, [$accOrder->[$i], $i] if($accOrder->[$i]);
	}
	my $numValidAccs = @validAccs;
	my $byteLen = ceil($numValidAccs * 2 / 8);
#	my $packFormat = sprintf "b%d", $numValidAccs * 2;
	my $packFormat = sprintf "b%d", $byteLen * 8;

	my($bin2bit, $bit2bin) = getBinBitMaps();

#	my $tmpDir = File::Temp->newdir(DIR => $self->{'tempdir'}, CLEANUP => 1);
	my $tmpDir = $self->{'tempdir'};

	my $packedFile = sprintf("%s/tasukePackedFile", $tmpDir);
	my $tmpBin;
	open($tmpBin, "> $packedFile");
	binmode $tmpBin;

#my @loop2chrpos;

	my $numSmplInfos = -1;
	my $loop = 0;
	my $in = _getFileHandle($vcf);
	while(<$in>) {
		chomp;
		my $datarow = $_;
		next if($datarow !~ m/\S/);
		next if($datarow =~ m/^#/);

		my @vals = split /\t/, $datarow;
		my $gtIdx = getVcfGtIndex($vals[8]);
		my @smplInfos = @vals[9 .. $#vals];
		my @alts = split /,/, $vals[4];
		my $bitStr = "";
		my $existsAfterFilter = 0;
		for(my $i=0; $i<$numValidAccs; $i++) {
			my $acc = $validAccs[$i]->[0];
			my $col = $validAccs[$i]->[1];

			my $bin;
			if($gtIdx < 0) {
				$bin = "NA";
			} else {
				my @smplVals = split /\:/, $smplInfos[$col];
				$bin = $self->getGtBinFilterVariantType($vals[3], \@alts, $smplVals[$gtIdx]);
				$existsAfterFilter = 1 if($bin eq '1' || $bin eq '2');
#				$bin = getGtBin($smplVals[$gtIdx]);
			}
			$bitStr .= $bin2bit->{$bin};
		}
#print STDERR "$bitStr\n";
#printf "%s\t%s\n", $vals[0], $vals[1];

		if($existsAfterFilter) {		
			print $tmpBin pack($packFormat, $bitStr);
			$loop++;
#push @loop2chrpos, [$vals[0], $vals[1]];

			printf STDERR "Progress: %s:%d (%d sites loaded)\n", $vals[0], $vals[1], $loop if(!($loop % 10000));		
		}
#last if($loop>=10000);
	}
	close($in);
	close($tmpBin);

	open($tmpBin, "<$packedFile");
	binmode($tmpBin);
	for(my $i=0; $i<$numValidAccs; $i++) {
		print $out $validAccs[$i]->[0];

#my $j = 0;
		sysseek($tmpBin, 0, 0);
		while(sysread $tmpBin, my $buf, $byteLen) {
			my $bit = '';
			$bit .= vec($buf, $i*2, 1);
			$bit .= vec($buf, $i*2 + 1, 1);
			printf $out ",%s", $bit2bin->{$bit};
#printf "%s\t%s\t%s\t%s\n", $loop2chrpos[$j]->[0], $loop2chrpos[$j]->[1], $validAccs[$i]->[0], $bit2bin->{$bit};
#$j++;
		}
		print $out "\n";

		printf STDERR "(%d/%d)\tOutput binary list about %s\t%d pos\n", ($i + 1), $numValidAccs, $validAccs[$i]->[0], $loop;
	}
	close($tmpBin);
}

###
 # For multi-sumple VCF file. All GT information extracted from the VCF is stored in memory.
###
sub outputFromMvcfMem {
	my $self = shift;
	my($vcf, $out, $accOrder) = @_;

	$accOrder = $self->getSamples($vcf) if(!$accOrder);	# Include null element.
	my @validAccs;
	my $numAccOrder = @{$accOrder};
	for(my $i=0; $i<$numAccOrder; $i++) {
		push @validAccs, [$accOrder->[$i], $i] if($accOrder->[$i]);
	}
	my $numValidAccs = @validAccs;

	my $packFormat = sprintf "b%d", $numValidAccs * 2;
	my($bin2bit, $bit2bin) = getBinBitMaps();

	my @acc2bins;
	for(my $i=0; $i<$numValidAccs; $i++) {
		push @acc2bins, [];	# Initialize
	}

	my @packedGtArr;
	my $loop = 0;
	my $in = _getFileHandle($vcf);
	while(<$in>) {
		chomp;
		my $datarow = $_;
		next if($datarow !~ m/\S/);
		next if($datarow =~ m/^#/);

		my @vals = split /\t/, $datarow;
		my $gtIdx = getVcfGtIndex($vals[8]);
		my @smplInfos = @vals[9 .. $#vals];
		my $bitStr = "";
		for(my $i=0; $i<$numValidAccs; $i++) {
			my $col = $validAccs[$i]->[1];
			my $bin;
			if($gtIdx < 0) {
				$bin = 'NA';
			} else {
				my @smplVals = split /\:/, $smplInfos[$col];
				$bin = getGtBin($smplVals[$gtIdx]);
			}
			$bitStr .= $bin2bit->{$bin};
		}
		$packedGtArr[$loop] = pack($packFormat, $bitStr);

		$loop++;
		if( !($loop % 10000) ) {
			printf STDERR "Progress: %s:%d (%d sites loaded)\n", $vals[0], $vals[1], $loop;
		}
#last if($loop >= 10000);
	}
	close($in);

	# Output
	print STDERR "Creating binary(0/1) table\n";
	for(my $i=0; $i<$numValidAccs; $i++) {
		my $acc = $validAccs[$i]->[0];
		print $out $acc;
		for(my $j=0; $j<$loop; $j++) {
			my $posBits = $packedGtArr[$j];

			my $bit = '';
			$bit .= vec($posBits, $i*2, 1);
			$bit .= vec($posBits, $i*2 + 1, 1);
			printf $out ",%s", $bit2bin->{$bit};
		}
		print $out "\n";

		printf STDERR "(%d/%d)\tOutput binary list about %s\t%d pos\n", ($i + 1), $numValidAccs, $acc, $loop;
	}
	undef @packedGtArr;
}

###
 # Methods for VCF.
 # A temporary DBM file is generated. It is slow, but can be used for data sets of around 1000 samples.
###
sub outputFromMvcfByDBM {
	my $self = shift;
	my($vcf, $out, $accOrder) = @_;

	$accOrder = $self->getSamples($vcf) if(!$accOrder);	# Include null element.
	my @validAccs;
	my $numAccOrder = @{$accOrder};
	for(my $i=0; $i<$numAccOrder; $i++) {
		push @validAccs, [$accOrder->[$i], $i] if($accOrder->[$i]);
	}
	my $numValidAccs = @validAccs;
	my $packFormat = sprintf "b%d", $numValidAccs * 2;

	my %gt2bit = (
		'0' => '00',
		'1' => '01',
		'2' => '10',
		'NA' => '11'
	);
	my %bit2gt;
	foreach my $gt (keys(%gt2bit)) {
		$bit2gt{$gt2bit{$gt}} = $gt;
	}

#	my $tmpDir = File::Temp->newdir(DIR => $self->{'tempdir'});
	my $tmpDir = $self->{'tempdir'};

	my $packedGtFile = sprintf("%s/packedGtArr", $tmpDir);
	my %packedGtHash;
	if($numValidAccs < 3000) {
		printf STDERR "Valid ACCs:%d\tUse SDBM_File\n", $numValidAccs;
		tie %packedGtHash, 'SDBM_File', $packedGtFile, O_RDWR|O_CREAT, 0666 or die "Couldn't tie SDBM file 'filename': $!; aborting";
	} else {
		printf STDERR "Valid ACCs:%d\tUse NDBM_File\n", $numValidAccs;
		tie %packedGtHash, 'NDBM_File', $packedGtFile, O_RDWR|O_CREAT, 0666 or die "Couldn't tie NDBM file 'filename': $!; aborting";
	}

	my $numSmplInfos = -1;
	my $in = _getFileHandle($vcf);
	my $loop = 0;
	while(<$in>) {
		chomp;
		my $datarow = $_;
		next if($datarow !~ m/\S/);

		if($datarow =~ m/^#/) {
			next;
		} else {
			my @vals = split /\t/, $datarow;
			my $gtIdx = getVcfGtIndex($vals[8]);
			my @smplInfos = @vals[9 .. $#vals];
			my $bitStr = "";
			for(my $i=0; $i<$numValidAccs; $i++) {
				my $acc = $validAccs[$i]->[0];
				my $col = $validAccs[$i]->[1];

				my $bin;
				if($gtIdx < 0) {
					$bin = "NA";
				} else {
					my @smplVals = split /\:/, $smplInfos[$col];
					$bin = getGtBin($smplVals[$gtIdx]);
				}
#print STDERR "$smplInfo -> $bin\n";
				$bitStr .= $gt2bit{$bin};
#print STDERR "$packFormat\t$bitStr\n";
			}
			$packedGtHash{$loop} = pack($packFormat, $bitStr);

			$loop++;
			if( !($loop % 10000) ) {
				printf STDERR "Progress: %s:%d (%d sites loaded)\n", $vals[0], $vals[1], $loop;
			}
		}
	}
	close($in);

	# Output
	print STDERR "Creating binary(0/1) table\n";
	for(my $i=0; $i<$numValidAccs; $i++) {
		my $acc = $validAccs[$i]->[0];
		print $out $acc;
		for(my $j=0; $j<$loop; $j++) {
			my $bitStr = unpack($packFormat, $packedGtHash{$j});
			my $gt = $bit2gt{substr($bitStr, $i * 2, 2)};
			print $out ",$gt";
		}
		print $out "\n";

		printf STDERR "(%d/%d)\tOutput binary list about %s\t%d pos\n", ($i + 1), $numValidAccs, $acc, $loop;
	}
}

###
 # Methods for VCF.
 # If you have a vcf index, no extra memory or temporary files are needed,
 # but it is very slow so cannot be used for large data set.
 # For several dozen samples.
###
sub outputFromBcftools {
	my $self = shift;
	my($vcf, $out, $bcftoolsPath, $accOrder) = @_;
	
	$accOrder = $self->getSamples($vcf) if(!$accOrder);	# Include null element.
	#TODO convert sample name.

	my $numAccOrder = @{$accOrder};
	for(my $i=0; $i<$numAccOrder; $i++) {
		next if(!$accOrder->[$i]);

		my $acc = $accOrder->[$i];
		print $out $acc;

		my $cmd = sprintf "%s view -O v -s %s %s | grep -v \"^#\" | cut -f 9,10", $bcftoolsPath, $acc, $vcf;
		my $loop = 0;
		open(CMD, "$cmd|");
		while(<CMD>) {
			chomp;
			my $datarow = $_;
			next if($datarow !~ m/\S/);

			my($format, $smplInfo) = split /\t/, $datarow;
			my $gtIdx = getVcfGtIndex($format);
			my $bin;
			if($gtIdx < 0) {
				$bin = "NA";
			} else {
				my @smplVals = split /\:/, $smplInfo;
				$bin = getGtBin($smplVals[$gtIdx]);
			}
			print $out ",$bin";
			$loop++;

			printf STDERR "Output binary table: %s %d sites\n", $acc, $loop if(!($loop % 10000));
		}
		close(CMD);

		print $out "\n";
	}
}

###
 # Standard methods for TASUKE.
 # Uses DBM and DB access. Although a temporary file is created, it uses little memory and is fast.
###
sub output {
	my $self = shift;
	my($out) = @_;

	my $lib = new CommonMysql();
	my $isMvcfMode = $lib->modeIsPlusMultiVcf($self->{'dbh'});

	$out = *STDOUT if(!$out);

	my $accRefs = $self->getAccessions();
	my $totalAccs = @{$accRefs};

	my $chrsSql = "";
	if(@{$self->{'chrs'}} > 0) {
		$chrsSql = "and chr in (";
		my @chrs = "";
		foreach my $chr (@{$self->{'chrs'}}) {
			$chrsSql .= "'$chr',";
		}
		$chrsSql =~ s/,\s*$//;
		$chrsSql .= ')';
	}

	#TODO It is correct to get the Variant and then compare and set it to "NA".
	my $filterSql = "";
#	$filterSql .= sprintf(" and QUAL>=%d", $self->{'min_qual'}) if($self->{'min_qual'});
#	$filterSql .= sprintf(" and DP>=%d", $self->{'min_dp'}) if($self->{'min_dp'});
#	$filterSql .= sprintf(" and (DP3+DP4)>=%d", $self->{'min_dp4'}) if($self->{'min_dp4'});

	# Make variant-type filter sql
	my $filterVtypeSql = "";
	if($self->{'targetv'}->{'s'} && $self->{'targetv'}->{'i'} && $self->{'targetv'}->{'d'}) {
		# Nothing to do.
	} else {
		$filterVtypeSql = " and variant in (";
		foreach my $vtype (keys(%{$self->{'targetv'}})) {
			$filterVtypeSql .= sprintf("'%s',", $vtype) if($self->{'targetv'}->{$vtype});
		}
		$filterVtypeSql =~ s/,$//;
		$filterVtypeSql .= ")";
	}
	print STDERR "filterVtypeSql: $filterVtypeSql\n";

	my $formatPack = 'L';	# L

	my %chr2posArray;
	my @chrs;
	{	
#		my $tmpDir = File::Temp->newdir(DIR => $self->{'tempdir'});
		my $tmpDir = $self->{'tempdir'};

		print STDERR "Start searching for Variant positions from DB\n";
#		my %uniAccChr;	# for Guide messages
		my %chr2pos;
		my $loop = 1;
		my $chrCnt = 0;
		foreach my $accRef (@{$accRefs}) {
			 printf STDERR "(%d/%d)\tScanning %s\n", $loop, $totalAccs, $accRef->{'accession_id'};
			
			my $sqlAcc;
			if($isMvcfMode) {
				$sqlAcc = sprintf "select chr,pos_vcf from %s_variant where GT!='0/0'", $accRef->{'accession_db'};
			} else {
				$sqlAcc = sprintf "select chr,IF(variant='d',pos-1,pos) as pos from %s_variant where GT!='0/0'", $accRef->{'accession_db'};
#				$sqlAcc = sprintf "select chr,pos from %s_variant where GT!='0/0'", $accRef->{'accession_db'};
			}
#			my $sqlAcc = sprintf "select chr,pos from %s_variant where GT!='0/0' and GT!='./.'", $accRef->{'accession_db'};

			$sqlAcc .= $chrsSql if($chrsSql);
			$sqlAcc .= $filterSql if($filterSql);
			$sqlAcc .= $filterVtypeSql if($filterVtypeSql);
#print STDERR "$sqlAcc\n";
			my $sthAcc = $self->{'dbh'}->prepare($sqlAcc);
			$sthAcc->execute();
			while(my $arrayref = $sthAcc->fetchrow_arrayref()) {
				my $chr = $arrayref->[0];
				my $pos = pack $formatPack, $arrayref->[1];	# To save memory, use the packed value(long-binary) as the key.
				if(!exists($chr2pos{$chr})) {
					my $dbfile = sprintf("%s/file%d", $tmpDir, $chrCnt);
					$chrCnt++;
					tie my %data, 'SDBM_File', $dbfile, O_RDWR | O_CREAT, 0666 or die "Couldn't tie SNDBM file 'filename': $!; aborting";
#tie my %data, 'GDBM_File', $dbfile, O_RDWR | O_CREAT, 0666 or die "Couldn't tie GNDBM file 'filename': $!; aborting";
#tie my %data, 'NDBM_File', $dbfile, O_RDWR | O_CREAT, 0666 or die "Couldn't tie NDBM file 'filename': $!; aborting";
#tie my %data, 'NDBM_File', $dbfile, O_RDWR | O_CREAT, 0666 or die "Couldn't tie NDBM file 'filename': $!; aborting";
					$chr2pos{$chr} = \%data;
				}
				$chr2pos{$chr}->{$pos} = undef;	# All you need is a key. undef is set to save memory.
			}
			$sthAcc->finish();

			$loop++
		}
#print STDERR "\%chr2pos:". total_size(\%chr2pos)."\n";

		@chrs = sort {$a cmp $b} keys(%chr2pos);
		if(1) {
			print STDERR "Sorting Chr+Pos\n";
			foreach my $chr (@chrs) {
				print STDERR "$chr\n";
				$chr2posArray{$chr} = [];

				my @poss;
				while ( my ( $key, $val ) = each %{$chr2pos{$chr}} ) {
					push @poss, int(unpack($formatPack, $key));
				}

				my @sorted = sort {$a <=> $b} @poss;
				$chr2posArray{$chr} = \@sorted;
				undef @poss;	# free up memory

				$chr2pos{$chr} = undef;	# free up memory
#print STDERR "\%chr2posArray:". total_size(\%chr2posArray). "   <- \%chr2pos:". total_size(\%chr2pos). "\n";
			}
		}
		undef %chr2pos;	# Since %chr2pos will not be used after that, free the memory.
	}

	if($self->{'check_depth_null'}) {
#		my $lib = new CommonMysql();
		print STDERR "Make Acc2Column2Table map\n";
		foreach my $accRef (@{$accRefs}) {
			$self->{'accdb2col2table'}->{$accRef->{'accession_db'}} = $lib->getCol2TableMap($self->{'dbh'}, $self->{'db'}, $accRef->{'accession_db'}. '_accession');
		}
	}

	# Output
	print STDERR "Creating binary(0/1) table\n";
	my $loop = 1;
	my $spliceUnit = 10000;
	my $progressBarResol = 30;
	foreach my $accRef (@{$accRefs}) {
		print $out $accRef->{'accession_id'};
		my $num = 0;
		foreach my $chr (@chrs) {
			my @VariantPoss = @{$chr2posArray{$chr}};	# Deep copy
			my $depthNullPoss = {};
			if($self->{'check_depth_null'} && $self->{'depth_null_in_memory'}) {
				my $colname = $chr. '_b';
				my $tablename = $self->{'accdb2col2table'}->{$accRef->{'accession_db'}}->{$colname};
				my @chrPos = @{$chr2posArray{$chr}};	# Deep copy
				my $unitCnt = ceil(@chrPos / $spliceUnit);
				my $loopChrAcc = 1;

				$| = 1;	# AutoFlush:On
				my $prevPogCount = 0;
				printf STDERR "\tGet Depth=NULL pos(%s)|", $chr;
				while(my @poss = splice(@chrPos, 0, $spliceUnit)) {
					my $progCount = ceil($loopChrAcc / $unitCnt * $progressBarResol);
					if($prevPogCount != $progCount) {
						print STDERR "*";
						$prevPogCount = $progCount;
					}
					my $whereSql = sprintf "pos in(%s)", join(',', @poss);
					my $sql = sprintf "select pos from %s where (`%s` is null or `%s`=0) and %s", $tablename, $colname, $colname, $whereSql;
					my $sthNull = $self->{'dbh'}->prepare($sql);
					$sthNull->execute();
					while(my $arrayref = $sthNull->fetchrow_arrayref()) {
						my $posPack = pack("L", $arrayref->[0]);
						$depthNullPoss->{$posPack} = undef;
					}
					$sthNull->finish();
					$loopChrAcc++;
				}
				print STDERR "\n";
				$| = 0;		# AutoFlush:Off
				undef @chrPos;
			}
#			my $depthNullPoss = $depthNullTable{$accRef->{'accession_db'}}->{$chr};

			my $sqlAcc;
			if($isMvcfMode) {
				$sqlAcc = sprintf "select pos_vcf,GT from %s_variant where GT!='0/0' and chr='%s'", $accRef->{'accession_db'}, $chr;
				$sqlAcc .= $filterSql if($filterSql);
				$sqlAcc .= $filterVtypeSql;
				$sqlAcc .= ' order by pos_vcf asc';
			} else {
				$sqlAcc = sprintf "select IF(variant='d',pos-1,pos) as pos,GT from %s_variant where GT!='0/0' and chr='%s'", $accRef->{'accession_db'}, $chr;
#				$sqlAcc = sprintf "select pos,GT from %s_variant where GT!='0/0' and chr='%s'", $accRef->{'accession_db'}, $chr;
				$sqlAcc .= $filterSql if($filterSql);
				$sqlAcc .= $filterVtypeSql;
				$sqlAcc .= ' order by pos asc';
			}
#			my $sqlAcc = sprintf "select chr,pos,GT from %s_variant where GT!='0/0' and chr='%s'", $accRef->{'accession_db'}, $chr;
#			$sqlAcc .= $filterSql if($filterSql);
#			$sqlAcc .= ' order by pos asc';

			my $sthAcc = $self->{'dbh'}->prepare($sqlAcc);
			$sthAcc->execute();
			while(my $arrayref = $sthAcc->fetchrow_arrayref()) {
				my $pos = $arrayref->[0];
				my $gt = $arrayref->[1];
				while(my $nextVariantPos = shift @VariantPoss) {	#TODO shift is costly.
					if($pos > $nextVariantPos) {
						# Pos where Variant was not detected
						if($self->{'check_depth_null'}) {
							if($self->{'depth_null_in_memory'}) {
								my $nextVariantPosPack = pack("L", $nextVariantPos);
								if(exists($depthNullPoss->{$nextVariantPosPack})) {
									print $out ",NA";
								} else {
									print $out ",0";	# 0/0:Ref Homo
								}
							} else {
								# Currently,  It will never come here.
								if($self->isPosNull($accRef->{'accession_db'}, $chr, $nextVariantPos)) {
									print $out ",NA";
								} else {
									print $out ",0";	# 0/0:Ref Homo
								}
							}
						} else {
							print $out ",0";	# 0/0:Ref Homo
#printf "%s\t%s\t%s\t0\n", $chr, $nextVariantPos, $accRef->{'accession_id'};
						}
						$num++;
					} else {
						# Pos where Variant was detected
						printf $out ",%s", getGtBin($gt);
#printf "%s\t%s\t%s\t%s\n", $chr, $nextVariantPos, $accRef->{'accession_id'}, getGtBin($gt);
						$num++;
						last;
					}
				}
			}
			$sthAcc->finish();
			while(my $nextVariantPos = shift @VariantPoss) {
				print $out ",0";
#printf "%s\t%s\t%s\t0\n", $chr, $nextVariantPos, $accRef->{'accession_id'};
				$num++;
			}
			undef $depthNullPoss;
		}
		printf STDERR "(%d/%d)\tOutput binary list about %s\t%d pos\n", $loop, $totalAccs, $accRef->{'accession_id'}, $num;
		$loop++;
		print $out "\n";
	}

	# Immediately release memory pool
	undef %chr2posArray;
}
sub isPosNull {
	my $self = shift;
	my($accdb, $chr, $pos) = @_;

	my $colname = sprintf "%s_b", $chr;
	my $tablename = $self->{'accdb2col2table'}->{$accdb}->{$colname};
	my $sql = sprintf "select `%s` from %s where pos=%d", $colname, $tablename, $pos;
	my $sth = $self->{'dbh'}->prepare($sql);
	$sth->execute();
	if(my $arrayref = $sth->fetchrow_arrayref()) {
		if($arrayref->[0]) {
			return 0;
		} else {
			# 0 or null
			return 1;
		}
	} else {
		return 0;
	}	
	$sth->finish();
}
sub getAccDbNum {
	my($accdb) = @_;
	if($accdb =~ m/^d(\d+)$/) {
		return $1;	
	} else {
		return 0;
	}
}
sub getAccessions {
	my $self = shift;

	my $accids = $self->{'accs'};

	my $accSql = "";
	if(@{$accids} > 0) {
		$accSql = 'accession_id in (';
		foreach my $accid (@{$accids}) {
			$accSql .= "'$accid',";
		}
		$accSql =~ s/,\s*$//;
		$accSql .= ')';
	}

	my @array;
	my $sql = "select * from sys_accession";
	$sql .= " where $accSql" if($accSql);
	my $sth = $self->{'dbh'}->prepare($sql);
	$sth->execute();
	while(my $hashref = $sth->fetchrow_hashref()) {
		push @array, $hashref;
	}
	$sth->finish();

	my @arraySorted  = sort {getAccDbNum($a->{'accession_db'}) <=> getAccDbNum($b->{'accession_db'})} @array;
	return \@arraySorted;
}
sub getChrs {
	my $self = shift;

	my $chrids = $self->{'chrs'};
	
	my $chrSql = "";
	if(@{$chrids} > 0) {
		$chrSql = 'id in (';
		foreach my $chrid (@{$chrids}) {
			$chrSql .= "'$chrid',";
		}
		$chrSql =~ s/,\s*$//;
		$chrSql .= ')';
	}

	my @array;
	my $sql = "select * from sys_chromosome order by seq";
	$sql .= " where $chrSql" if($chrSql);
	my $sth = $self->{'dbh'}->prepare($sql);
	$sth->execute();
	while(my $hashref = $sth->fetchrow_hashref()) {
		push @array, $hashref;
	}
	$sth->finish();

	return \@array;	
}
sub getChrsMap {
	my $self = shift;

	my $chrs = $self->getChrs();
	my %map;
	foreach my $hashref (@{$chrs}) {
		$map{$hashref->{'id'}} = $hashref;
	}
	return \%map;
}
sub getChrIds {
	my $self = shift;

	my $chrRefs = $self->getChrs();
	my @array;
	foreach my $chrRef (@{$chrRefs}) {
		push @array, $chrRef->{'id'};
	}
	return \@array;
}
sub getLongestChrId {
	my $self = shift;

	my $chrRefs = $self->getChrs();
	my $longestChrid = "";
	my $maxLength = 0;
	foreach my $chrRef (@{$chrRefs}) {
		if(!$longestChrid || $maxLength < $chrRef->{'size'}) {
			$longestChrid = $chrRef->{'id'};
			$maxLength = $chrRef->{'size'};
		}
	}
	return $longestChrid;
}
sub makeNewickByPhylip {
    my $self = shift;
    my($phylipDirPath, $matrixPath, $newickPath) = @_;

    my $currentDir = getcwd();

    my $stat = 1;
    eval {
#	my $wkdir = File::Temp->newdir(DIR => $self->{'tempdir'}, CLEANUP => 1);
	my $wkdir = $self->{'tempdir'}. '/phylip';
	mkdir $wkdir or die "Temp-dir for phylip failed: $wkdir\n";

	my $cmd;

	# exec Phylip neighbor.
	$cmd = sprintf "ln -s %s %s/infile", File::Spec->rel2abs($matrixPath), $wkdir;
	printf STDERR "CMD: %s\n", $cmd;
	system($cmd) == 0 or die "Can't create phylip neighbor infile: $!\n";

	chdir($wkdir);
	$cmd = sprintf("echo -e \"L\nY\n\" | %s/exe/neighbor", $phylipDirPath);
	printf STDERR "CMD: %s\n", $cmd;
	system($cmd) == 0 or die "Phylip neighbor failed: $!\n";

	# Mid-point rooting
	my $phylipRetreePath = sprintf "%s/exe/retree", $phylipDirPath;
	if(-f $phylipRetreePath) {
		$cmd = sprintf("echo -e \"Y\nouttree\nM\nW\nR\nR\nQ\n\" | %s > /dev/null", $phylipRetreePath);
		printf STDERR "CMD: %s\n", $cmd;
		system($cmd) == 0 or die "Phylip retree failed: $!\n";
	} else {
		printf STDERR "Phylip retree(%s) not found. No Mid-point rooting will be performed.\n", $phylipRetreePath;
	}

	chdir($currentDir);

	$cmd = sprintf("mv -f %s/outtree %s", $wkdir, $newickPath);
	printf STDERR "CMD: %s\n", $cmd;
	system($cmd) == 0 or die "Can't move newick file: $!\n";
	$stat = 0;
    };
    if ($@) {
	print STDERR "Exception occur: $@";
	chdir($currentDir);
    }
    return $stat;
}
sub makeNewickByPhylipiViaRename {
    my $self = shift;
    my($phylipDirPath, $matrixPath, $newickPath) = @_;

    my $stat = 1;
    eval {

	my $afNamePat = "x%09d";
	my $smplIdx = 0;
	my $encountSmpl = 1;

	my($AF, $afFile) = tempfile(DIR => $self->{'tempdir'}, UNLINK => 1);
	my $BF;
	open($BF, "<$matrixPath") or die "Distance matrix open error: $!";
	sysseek $BF, 0, 0;

	# Print header(sumple total)
	while(sysread $BF, my $bin, 1) {
		my $char = unpack('a1', $bin);
		print $AF $char;
		last if($char eq "\n");
	}

	# Print body while changing the sample name.
	my %tmp2ori;
	while(sysread $BF, my $bin, 1) {
		my $char;
		if($encountSmpl) {
			my $tmpName = sprintf $afNamePat, $smplIdx;
			print $AF $tmpName;
			$smplIdx++;
			my $oriName = unpack('a1', $bin);
			while(sysread $BF, my $binsub, 1) {
				$char = unpack('a1', $binsub);
				if($char eq " " || $char eq "\n") {
					$tmp2ori{$tmpName} = $oriName;
					last;
				}
				$oriName .= $char;
			}
			printf STDERR "INFO: Replace ID in distance matrix file: %s -> %s\n", $oriName, $tmpName;
			$encountSmpl = 0;
		} else {
			$char = unpack('a1', $bin);
		}
		print $AF $char;

		$encountSmpl = 1 if($char eq "\n");
	}
	close($BF);
	close($AF);

	my ($preNewickPathFh, $preNewickPath) = tempfile(DIR => $self->{'tempdir'}, UNLINK => 1); 
	close($preNewickPathFh);
	my $stat = $self->makeNewickByPhylip($phylipDirPath, $afFile, $preNewickPath);

	# Replace temp-ID in Newick with the real-ID. Remove line breaks(\n).
	my $newickStr;
	open($BF, "<$preNewickPath") or die "Pre-newick file open error: $!";
	while(<$BF>) {
		chomp;
		my $datarow = $_;
		next if($datarow !~ m/\S/);

		$newickStr .= $datarow;
	}
	close($BF);
	foreach my $tmpName (keys(%tmp2ori)) {
		my $repStat = $newickStr =~ s/$tmpName/$tmp2ori{$tmpName}/;
		printf STDERR "INFO: RestoreID in  newick file: %s -> %s = %s\n", $tmpName, $tmp2ori{$tmpName}, $repStat;
	}
	open($AF, ">$newickPath") or die "Finally newick file open error: $!";
	print $AF $newickStr;
	close($AF);

	$stat = 0;
    };
    if ($@) {
    	print STDERR "Exception occur: $@";
    }

    return $stat;
}

1;
