package ChemicalReaction;
use Common;
@ISA = qw(Common);

use strict;
#use warnings;

use Math::Matrix;
use Math::MatrixReal;

use CSV;
use Crystal::AtomType;

#===============================================
# コンストラクタ、デストラクタ
#===============================================
sub new
{
	my ($module, $pDBFiles, $BaseDBDir) = @_;
	my $this = {};
	bless $this;

	$this->Initialize();

	$this->ReadDBFiles($pDBFiles, $BaseDBDir) if($pDBFiles ne '');
	return $this;
}

sub DESTROY
{
	my ($this) = @_;
	$this->SUPER::DESTROY(@_);
}

sub Initialize
{
	my ($this, $InitializeDB) = @_;
	$InitializeDB = 1 if(!defined $InitializeDB);

	$this->ClearHashes();
	if($InitializeDB) {
		delete $this->{pDBs};
	}
}

sub ClearHashes
{
	my ($this) = @_;
	$this->{pReagents}    = [];
	$this->{Reagents}     = [];
	$this->{ReagentsHash} = {};
	$this->{pProducts}    = [];
	$this->{Products}     = [];
	$this->{ProductsHash} = {};
}

#===============================================
# 一般関数
#===============================================
sub pDBs            { my ($this)=@_; return $this->{pDBs}; };
sub DB              { my ($this, $index)=@_; return $this->{pDBs}->[$index]; };
sub InitialReaction { my ($this)=@_; return $this->{Initial} ; };
sub FinalReaction   { my ($this)=@_; return $this->{Final} ; };
sub Reagents        { my ($this)=@_; return $this->{Initial} ; };
sub Products        { my ($this)=@_; return $this->{Final} ; };
sub ReagentsHash    { my ($this)=@_; return $this->{ReagentsHash} ; };
sub ProductsHash    { my ($this)=@_; return $this->{ProductsHash} ; };
sub pCompounds {
	my ($this, $key) = @_;
	if($key =~ /^P/i) {
		return $this->{ProductsHash}->{pCompounds};
	}
	return $this->{ReagentsHash}->{pCompounds};
}
sub nCompound {
	my ($this, $key, $compound) = @_;
	
	my $pHash;
	if($key =~ /^P/i) {
		$pHash = $this->ProductsHash();
	}
	else {
		$pHash = $this->ReagentsHash();
	}
	return $pHash->{$compound}{n};
}
sub ElementsArray {
	my ($this, $key) = @_;
	return $this->pElements($key);
}
sub pElements {
	my ($this, $key) = @_;
	if($key =~ /^P/i) {
		return $this->{ProductsHash}->{pElements};
	}
	return $this->{ReagentsHash}->{pElements};
}
sub nElement {
	my ($this, $key, $element) = @_;

	my $pHash;
	if($key =~ /^P/i) {
		$pHash = $this->ProductsHash();
	}
	else {
		$pHash = $this->ReagentsHash();
	}
	return $pHash->{$element}{nElement};
}
sub MolcularWeight
{
	my ($this, $Compound) = @_;

	my $R = new ChemicalReaction;
	$R->AnalyzeAReaction('Reagents', $Compound);
	my $pHash = $R->ReagentsHash();
	my $pElements = $R->ElementsArray("Reagents");
	my $w = 0.0;
	for(my $i = 0 ; $i < @$pElements ; $i++) {
		my $e = $pElements->[$i];
		my $n = $R->nElement("Reagents", $e);
#print "$e: $n\n";
		my $pAtomDB = AtomType::GetAtomDBHash($e);
#print "  Z=$pAtomDB->{Z}: M=$pAtomDB->{mass}\n";
		$w += $n * $pAtomDB->{mass};
	}
	return $w;
}

sub IsIncluded
{
	my ($this, $key, $pArray, $CaseSensitive) = @_;
	$CaseSensitive = 0 if(!defined $CaseSensitive);

	$key = uc $key if(!$CaseSensitive);
	for(my $i = 0 ; $i < @$pArray ; $i++) {
		if($CaseSensitive) {
			return 1 if($pArray->[$i] eq $key);
		}
		else {
			return 1 if(uc $pArray->[$i] eq $key);
		}
	}
	return 0;
}

sub SplitReaction
{
	my ($this, $Reaction) = @_;
	return Utils::Split("=>?", $Reaction);
}

sub ReactionToCompounds
{
	my ($this, $reaction) = @_;

	my @Chemicals = Utils::Split("[\\+]", $reaction);
	my (@Compounds, @n);
	for(my $i = 0 ; $i < @Chemicals ; $i++) {
		($n[$i], $Compounds[$i]) = ($Chemicals[$i] =~ /^([\d\.\/]*)(\D.*)$/);
		$n[$i] = 1 if($n[$i] eq '');
		$n[$i] = eval($n[$i]);
	}
	return (\@Chemicals, \@Compounds, \@n);
}

sub CompoundToElements
{
	my ($this, $Compound) = @_;

	$Compound =~ s/:.*$//;
	$Compound =~ s/\-.*$//;
	$Compound =~ s/\(.*(AFM|FM|NM|NP|SP|NM|atom).*\)$//i;
	my (@a) = Utils::FindAll("[A-Z][a-z]*\\d*", $Compound);
	my (@Element, @n);
	for(my $i = 0 ; $i < @a ; $i++) {
		($Element[$i], $n[$i]) = ($a[$i] =~ /^([A-Z][a-z]*)(\d*)$/);
		$n[$i] = 1 if(!defined $n[$i] or $n[$i] eq '');
	}
	return (\@Element, \@n);
}

sub AnalyzeAReaction
{
	my ($this, $key, $reaction) = @_;

	my ($pChemicals, $pCompound, $pnCompound) = $this->ReactionToCompounds($reaction);

	my ($pReaction, $pList, $pHash);
	if($key =~ /^P/i or $key =~ /^F/i) {
		$this->{pProducts} = $pChemicals;
		$pList     = $this->{Products};
		$pHash     = $this->{ProductsHash};
	}
	else {
		$this->{pReagents} = $pChemicals;
		$pList     = $this->{Reagents};
		$pHash     = $this->{ReagentsHash};
	}

#print "[$key]:\n";
	my %ElementsHash;
	for(my $i = 0 ; $i < @$pCompound ; $i++) {
#print "R[$i]=$pChemicals->[$i]\n";
		my $compound = $pCompound->[$i];
		my $n        = $pnCompound->[$i];
		$pList->[$i]{n}        = $n;
		$pList->[$i]{Compound} = $compound;
		$pHash->{$compound}{n} = $n;
#print "  ${n}[$compound]: ";

		my ($pElement, $pnElement) = $this->CompoundToElements($compound);
		for(my $i = 0 ; $i < @$pElement ; $i++) {
			my $element = $pElement->[$i];
			my $ne      = $pnElement->[$i];
#print "  $ne * $element + ";
			$ElementsHash{$element}++;
			$pHash->{$element}{ne} = $ne;
			$pHash->{$compound}{$element}{ne} += $ne;
			$pHash->{$element}{nElement} += $n * $ne;
			$pHash->{nElement}{$element} += $n * $ne;
		}
#print "\n";
	}
#print "E: ", join(', ', keys %ElementsHash), "\n";
	$pHash->{pCompounds} = $pCompound; #\@Compounds;
	$pHash->{pElements} = [keys %ElementsHash];
}

sub Analyze
{
	my ($this, $Reaction, $IsPrint) = @_;
	$IsPrint = 1 if(!defined $IsPrint);

	if($Reaction eq 'list') {
		$this->ListCompounds();
		return ();
	}
	elsif($Reaction eq 'exit' or $Reaction eq 'quit' ) {
		exit;
	}

	$this->{Reaction} = $Reaction;
	($this->{Initial}, $this->{Final}) = $this->SplitReaction($Reaction);

#	$this->{Initial} =~ s/\-/+\-/g;
#	$this->{Final}   =~ s/\-/+\-/g;
	$this->AnalyzeAReaction("Reagents", $this->{Initial});
	$this->PrintElements("Reagents") if($IsPrint);
	$this->AnalyzeAReaction("Products", $this->{Final});
	$this->PrintElements("Products") if($IsPrint);
}

sub ListCompounds
{
	my ($this) = @_;
	my @a = $this->SearchDBHit(".*");
	for(my $i = 0 ; $i < @a ; $i++) {
		print "$i: $a[$i]\n";
	}
}

sub CheckMassBalance
{
	my ($this, $IsPrint) = @_;
	$IsPrint = 0 if(!defined $IsPrint);

	$this->PrintnElements("Reagents") if($IsPrint);
	$this->PrintnElements("Products") if($IsPrint);

	my $pnInitialElementHash = $this->{ReagentsHash}{nElement};
	my $pnFinalElementHash   = $this->{ProductsHash}{nElement};
	foreach my $e (sort keys %$pnInitialElementHash) {
		if($pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}) {
print("Error: Mass balance violated for [$e]: "
	."$pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}\n") if($IsPrint);
			return 0;
		}
	}
	foreach my $e (sort keys %$pnFinalElementHash) {
		if($pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}) {
print("Error: Mass balance violated for [$e]: "
	."$pnInitialElementHash->{$e} != $pnFinalElementHash->{$e}\n") if($IsPrint);
			return 0;
		}
	}
	return 1;
}

sub GetDBHash
{
	my ($this, $Compound) = @_;

	for(my $i = 0 ; ; $i++) {
		my $DB = $this->DB($i);
		last if(!$DB);

		my $pHash = $DB->HashByXDataKey($Compound);
		return $pHash if($pHash and $pHash->{Status} eq '');
	}
	return undef;
}

sub SearchDBHit
{
	my ($this, $RegExp) = @_;

	my @Hit;
	for(my $i = 0 ; ; $i++) {
		my $DB = $this->DB($i);
		last if(!$DB);

		my $pX = $DB->{pX};
		$pX = $DB->GetData(0) if(!defined $pX);
		my $pS = $DB->GetData("Status");
#print "ps=$pS\n";
		my $nData = @$pX;
		for(my $i = 0 ; $i < $nData ; $i++) {
			if($pX->[$i] =~ /$RegExp/) {
				next if($pS and $pS->[$i] ne '');
				push(@Hit, $pX->[$i]);
			}
		}
	}
	return @Hit;
}

sub CalEForm
{
	my ($this) = @_;

	my @Reagents = @{$this->{pReagents}};
	my @Products = @{$this->{pProducts}};
	my $csv = $this->DB();

	my $CouldNotFind = 0;
	my $EFormInitial = 0;
	for(my $i = 0 ; $i < @Reagents ; $i++) {
		my $n = $this->{Reagents}[$i]{n};
		my $c = $this->{Reagents}[$i]{Compound};
		my $pHash = $this->GetDBHash($c);
		if(!defined $pHash) {
print "Error: Can not find [$c] in DBs.\n";
			$this->ShowSpeculatedSpecies($c);
			$CouldNotFind = 1;
			next;
#			return ();
		}
		my $Z    = $pHash->{Z};
		my $Etot = $pHash->{Etot};
		if(!defined $Etot) {
print "Error: Can not find [$c] in DBs.\n";
			$this->ShowSpeculatedSpecies($c);
			$CouldNotFind = 1;
			next;
#			return ();
		}
		$EFormInitial += $n * $Etot / $Z;
	}

	my $EFormFinal = 0;
	for(my $i = 0 ; $i < @Products ; $i++) {
		my $n = $this->{Products}[$i]{n};
		my $c = $this->{Products}[$i]{Compound};
		my $pHash = $this->GetDBHash($c);
		if(!defined $pHash) {
print "Error: Can not find [$c] in DBs.\n";
			$this->ShowSpeculatedSpecies($c);
			$CouldNotFind = 1;
			next;
#			return ();
		}
		my $Z    = $pHash->{Z};
		my $Etot = $pHash->{Etot};
		if(!defined $Etot) {
print "Error: Can not find [$c] in DBs.\n";
			$this->ShowSpeculatedSpecies($c);
			$CouldNotFind = 1;
			next;
#			return ();
		}
		$EFormFinal += $n * $Etot / $Z;
	}
	if($CouldNotFind) {
		return ();
	}
	else {
		my $EForm = $EFormFinal - $EFormInitial;
#print "EForm         : $EForm eV ($EFormFinal - $EFormInitial)\n";
		return ($EForm, $EFormInitial, $EFormFinal);
	}
}

sub ShowSpeculatedSpecies
{
	my ($this, $Compound) = @_;
	my @Hit = $this->SearchDBHit("^${Compound}[\\-:]");
	for(my $i = 0 ; $i < @Hit ; $i++) {
		print "   $Hit[$i]\n";
	}
}

sub ReadDBFile
{
	my ($this, $DBFile, $BaseDBDir) = @_;
	$this->{BaseDBDir} = $BaseDBDir if(defined $BaseDBDir);
print "ReadDBFile [$DBFile]\n";

	my $csv = new CSV;
	$DBFile = Deps::MakePath($BaseDBDir, $DBFile, 0) if(defined $BaseDBDir);
	if(!$csv->Read($DBFile, 0, 1)) {
print "Error: $!: Can not read [$DBFile].\n";
		return undef;
	}
	my $pX = $csv->GetXData("Compound");
	if(!$pX) {
		print "Error: Can not find 'Compound' tag in [$DBFile].\n";
	}
	for(my $i = 0 ; $i < @$pX ; $i++) {
		$pX->[$i] =~ s/Atom/atom/gi;
		$pX->[$i] =~ s/\((AFM|FM|NM|NP|SP|NM|atom)\)/:$1/gi;
		$pX->[$i] =~ s/([A-Z].*)-/$1:/g;
	}

	return $csv;
}

sub ReadDBFiles
{
	my ($this, $pDBFiles, $BaseDBDir) = @_;

	my @DBs;
	for(my $i = 0 ; $i < @$pDBFiles ; $i++) {
#print "f=$pDBFiles->[$i]\n";
		$DBs[$i] = $this->ReadDBFile($pDBFiles->[$i], $BaseDBDir);
		if(!$DBs[$i]) {
#			return undef;
		}
		else {
print "Success: DB $pDBFiles->[$i] has been loaded.\n";
		}
	}
	$this->{pDBs} = \@DBs;
	return $this->{pDBs};
}

sub PrintElements
{
	my ($this, $key) = @_;

	my ($pReaction, $pList, $pHash);
	if($key =~ /^P/i or $key =~ /^F/i) {
		$pReaction = $this->{pProducts};
		$pList     = $this->{Products};
		$pHash     = $this->{ProductsHash};
	}
	else {
		$pReaction = $this->{pReagents};
		$pList     = $this->{Reagents};
		$pHash     = $this->{ReagentsHash};
	}

	print "$key:\n";
	for(my $i = 0 ; $i < @$pReaction ; $i++) {
		my $n        = $pList->[$i]{n};
		my $compound = $pList->[$i]{Compound};
		print "  $n * $compound: ";

		my $IsFirst = 1;
		my $pElements = $pHash->{pElements};
		foreach my $element (sort @$pElements) {
#			my $ne        = $pHash->{$element}{ne};
			my $ne        = $pHash->{$compound}{$element}{ne};;
			next if($ne eq '');
			my $nElement  = $pHash->{$element}{nElement};

			print " + " if(!$IsFirst);
			print "$ne * $element";
			$IsFirst = 0;
		}
		print "\n";
	}

}

sub PrintnElements
{
	my ($this, $key) = @_;

	my $pnElementHash;
	if($key =~ /^P/i or $key =~ /^F/i) {
		$pnElementHash = $this->{ProductsHash}{nElement};
	}
	else {
		$pnElementHash = $this->{ReagentsHash}{nElement};
	}

	my $IsFirst = 1;
	print "nElements in $key: ";
	foreach my $e (sort keys %$pnElementHash) {
		print " + " if(!$IsFirst);
		print "$pnElementHash->{$e} * $e";
		$IsFirst = 0;
	}
	print "\n";
}
1;