#
# 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 Data::Dumper;
use POSIX qw(ceil);
use CommonLibTASUKE;

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,
		@_
	};
	return bless $self, $class;
}
sub verbose {
	my $self = shift;
	my($flg) = @_;
	$self->{'verbose'} = $flg;
}
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 output {
	my $self = shift;
	my($out) = @_;

	$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'});

	print STDERR "Start searching for Variant positions from DB\n";
#	my %uniAccChr;	# for Guide messages
	my %chr2pos;
	my $loop = 1;
	foreach my $accRef (@{$accRefs}) {
		 printf STDERR "(%d/%d)\tScanning %s\n", $loop, $totalAccs, $accRef->{'accession_id'};
		my $sqlAcc = sprintf "select chr,pos from %s_variant where GT != '0/0'", $accRef->{'accession_db'};
		$sqlAcc .= $chrsSql if($chrsSql);
		$sqlAcc .= $filterSql if($filterSql);
#print STDERR "$sqlAcc\n";
		my $sthAcc = $self->{'dbh'}->prepare($sqlAcc);
		$sthAcc->execute();
		while(my $arrayref = $sthAcc->fetchrow_arrayref()) {
			my $chr = $arrayref->[0];
			my $pos = $arrayref->[1];
			$chr2pos{$chr}->{$pos} = undef;	# All you need is a key. undef is set to save memory.
		}
		$sthAcc->finish();

		$loop++
	}

	my %chr2posArray;
	my @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 = keys(%{$chr2pos{$chr}});
			my @sorted = sort {$a <=> $b} @poss;
			$chr2posArray{$chr} = \@sorted;

			# free up memory
			undef @poss;
			$chr2pos{$chr} = undef;
		}
	}
	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";
	$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 $progressBar = '*' x ceil($loopChrAcc / $unitCnt * $progressBarResol);
#					printf STDERR "\tGet Depth=NULL pos(%s)|%-". $progressBarResol. "s|", $chr, $progressBar;	#Return carriage to the left end by \r.
					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()) {
						$depthNullPoss->{$arrayref->[0]} = undef;
					}
					$sthNull->finish();
					$loopChrAcc++;
#					print STDERR "\r";
#my $t = 0;
#for(my $i=0; $i<50000; $i++) { $t += $i;}
				}
				print STDERR "\n";
				$| = 0;		# AutoFlush:Off
				undef @chrPos;
			}
#			my $depthNullPoss = $depthNullTable{$accRef->{'accession_db'}}->{$chr};

			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()) {
				while(my $nextVariantPos = shift @VariantPoss) {	#TODO shift is costly.
					my $chr = $arrayref->[0];
					my $pos = $arrayref->[1];
					my $gt = $arrayref->[2];
					if($pos > $nextVariantPos) {
						# Pos where Variant was not detected
						if($self->{'check_depth_null'}) {
							if($self->{'depth_null_in_memory'}) {
								if(exists($depthNullPoss->{$nextVariantPos})) {
									print $out ",NA";
								} else {
									print $out ",0";	# 0/0:Ref Homo
								}
							} else {
								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
						}
						$num++;
					} else {
						# Pos where Variant was detected
						my @gtParam = split(/[\/\|]/, $gt);
						if($gtParam[0] eq '0' || (@gtParam > 1 && $gtParam[1] eq '0')) {
							print $out ",2";	# 0/1:Hetero
						} else {
							print $out ",1";	# 1/1:Alt Homo
						}
						$num++;
						last;
					}
				}
			}
			$sthAcc->finish();
			while(my $nextVariantPos = shift @VariantPoss) {
				print $out ",0";
				$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 %depthNullTable;
	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;
}

1;
