#!/usr/bin/perl

BEGIN {
use lib 'd:/Programs/Perl/lib';
use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;

use MyApplication;
use JFile;
use Sci::ChemicalReaction;
use Sci::Optimize;

#my $BaseDBDir = "d:/Programs.old/Perl/ChemicalReaction/Databases";
my $BaseDBDir = "./Databases";
#my $LogFile   = "d:/Programs/Perl/ChemicalReaction/Reaction-out.txt";
my $LogFile   = "./Reaction-out.txt";
#my $CheckFile = "d:/Programs/Perl/ChemicalReaction/Reaction-check.txt";
my $CheckFile = "./Reaction-check.txt";
my @DBFiles   = ("TotalEnergy.csv", 
#			"AeFe2As2-High.csv",
#			"AeCuFCh-High.csv", "AeCuFCh.csv", 
#			"LaMnOPn.csv", 
#			"BaFe2As2.csv", "SrFe2As2.csv", 
#			"AeFe2As2-High.csv",
#			"LaCuOSe.csv", 
#			"LaZnOP-all.csv",
#			"AB-all.csv",
#			"aIGZO-Dens608.csv", 
#			"ZnO.csv",
#			"aIGZO-Dens608-all.csv", "AOS-PBE96-all.csv", "C12A7-all.csv",
#			"LaMnOPn-all.csv", "LaTMOPn-AFM-Collinear-new2-all.csv",
#			"LaTMOPn-AFM-Collinear-new-all.csv", "LaTMOPn-AFM-Collinear-all.csv",
#			"Root-all.csv",
			);
my @RemoveCompounds = ();
#my @RemoveCompounds = ("FeO", "LaAs");

#my $Method = "PDL::Simplex";
my $Method = "Amoeba::Simplex";
#my $Method = "ModifiedNewton";
#my $Method = "LinearOptimization"; # 線形パラメータのみ。

my $EPS         = 1.0e-8;
my $nMaxIter    = 100000;
my $iPrintLevel = 1;
my $iPrintInterval = 20;
my $iIteration  = 0;
my $KScale0     = 0.1;
my $KScale      = 1.0;

print "use $Method\n";

my $PrintDetail = 0;
my $Reaction    = ($ARGV[0] eq '')? "La+Fe+O+As-0.1As" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "Sr+2Fe+2As+0.01O+0.01H" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "0.9Eu+Cu+F+Se" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "0.2H2O+0.1BaO" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "Ba+2Fe+2As+0.1H2" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "2Cu+F2" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "Se+Cu+2La+Mn+P+O2" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "SrF2+SrSe+Cu2Se+La+Mn+P+O2" : $ARGV[0];
#my $Reaction    = ($ARGV[0] eq '')? "CaF2+2CaS+2Cu2S" : $ARGV[0];

my $log = new JFile;
$log->Open($LogFile, "a") or die "$!: Can not append to [$LogFile].\n";
my $check = new JFile;
$check->Open($CheckFile, "w") or die "$!: Can not write to [$CheckFile].\n";

my $Opt = new Optimize;

#=============================================
# Create ChemicalReaction object
#=============================================
my $R = new ChemicalReaction(\@DBFiles, $BaseDBDir, "Composition");
if(!$R) {
	print "Error: $!: Can not create ChemicalReaction Object.\n";
	exit;
}
if(!$R->GetCSVDB()) {
	print "Error: $!: Can not read a DB.\n";
	exit;
}
$R->SetMode('Composition');
$R->PrintDBData(StopByError => 1);

#my $pCompoundDB = $R->GetCompoundDB();
my $pCompoundDB = $R->GetCompositionCompoundDB();

#=============================================
# Analyze reaction
#=============================================
print "Reaction: $Reaction\n";
$R->Analyze($Reaction, 0);
my $pReagentHash      = $R->ReagentsHash();
my $pReagentCompounds = $R->pCompounds('Reagents');
my $pReagentElements  = $R->pElements("Reagents");
# Reagent中の組成
my %nReagentElement;
for(my $i = 0 ; $i < @$pReagentElements ; $i++) {
	my $e = $pReagentElements->[$i];
	my $n = $R->nElement("Reagents", $e);
	$nReagentElement{$e} = $n;
	my $pAtomDB = AtomType::GetAtomDBHash($e);
	&mprint($check, "$e: n=$n: Z=$pAtomDB->{Z}: M=$pAtomDB->{mass}\n");
}

my $pElements    = $R->pElements("Reagents");
my @PossibleCompounds = $R->ExtractPossibleCompounds($pElements, [keys %$pCompoundDB], $pCompoundDB);
@PossibleCompounds = $R->RemoveCompounds(\@PossibleCompounds, \@RemoveCompounds);
if(@RemoveCompounds > 0) {
	$log->print("Removed: ", join(',', @RemoveCompounds), "\n");
	&mprint($check, "Removed: ", join(',', @RemoveCompounds), "\n");
}
&mprint($check, "To be considered: Atom\n");
foreach my $c (sort @PossibleCompounds) {
	next if($pCompoundDB->{$c}->{Type} ne 'Atom');
	&mprintf($check, "  %-2s %6.3f eV", $c, $pCompoundDB->{$c}->{Emol});
	&mprint($check, " ($pCompoundDB->{$c}->{Compound} in $pCompoundDB->{$c}->{Ref})\n");
}
print "To be considered: Elementary\n";
foreach my $c (sort @PossibleCompounds) {
	next if($pCompoundDB->{$c}->{Type} ne 'Elementary');
	&mprintf($check, "  %-2s %6.3f eV", $c, $pCompoundDB->{$c}->{Emol});
	&mprint($check, " $pCompoundDB->{$c}->{Type} ($pCompoundDB->{$c}->{Compound} in $pCompoundDB->{$c}->{Ref})\n");
}
&mprint($check, "To be considered: Compound\n");
foreach my $c (sort @PossibleCompounds) {
	next if($pCompoundDB->{$c}->{Type} ne 'Compound');
	&mprintf($check, "  %s %10.3f eV", $c, $pCompoundDB->{$c}->{Emol});
	&mprint($check, " $pCompoundDB->{$c}->{Type} ($pCompoundDB->{$c}->{Compound} in $pCompoundDB->{$c}->{Ref})\n");
}

my $EReaction = $R->CalTotalEnergy($Reaction);
&mprint($check, "Reaction: $Reaction\n");
&mprint($check, "  E(Reaction)=$EReaction eV\n");
$log->print("Reaction: $Reaction\n");
$log->print("  E(Reaction)=$EReaction eV\n");

my ($min, $MaxNInReaction) = Utils::CalMinMax([values %nReagentElement]);
&mprint($check, "\n");
&mprint($check, "Max nElement for [$Reaction]: $MaxNInReaction\n");

print "\n";
my @Results;
my $c = 0;
for(my $l = -2 ; $l < @PossibleCompounds ; $l++) {
#for(my $l = -2 ; $l <= -2 ; $l++) {
print "Optimization iter l=$l: MaxNInReaction=$MaxNInReaction\n";
	my (@Guess, @OptId, @Scale);
#	srand(1.0);

	my $pComposition;
	my $k = $KScale0;
	if(0 <= $l and $l < @PossibleCompounds) {
		for(my $i = 0 ; $i < @Guess ; $i++) {
			$pComposition->[$i] = ($i == $l)? $MaxNInReaction : 0.0;
		}
		$pComposition = $R->RepairComposition($pComposition, $pReagentElements, \@PossibleCompounds, $pCompoundDB, $PrintDetail);
		$k = $KScale0;
	}
	elsif($l == -1) {
		for(my $i = 0 ; $i < @Guess ; $i++) {
			$pComposition->[$i] = $MaxNInReaction;
		}
		$pComposition = $R->RepairComposition($pComposition, $pReagentElements, \@PossibleCompounds, $pCompoundDB, $PrintDetail);
		$k = $KScale0;
	}
	else {
		$pComposition = $R->MakePossibleComposition($pReagentElements, \@PossibleCompounds, $pCompoundDB, $MaxNInReaction, $PrintDetail);
		$k = $KScale;
	}
	@Guess = @$pComposition;
	for(my $i = 0 ; $i < @Guess ; $i++) {
		$OptId[$i] = 1;
		$Scale[$i] = $k * $MaxNInReaction;
	}

#必要なければ、PDLDisplayFunc以降は指定する必要はない
	my ($pMinComposition, $MinE) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				\&MinimizationFunc, \&PDLDisplayFunc,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
#print "Optimized at (",join(',',@{$pMinComposition}),")=$MinEform\n";
	$R->RepairComposition($pMinComposition, $pReagentElements, \@PossibleCompounds, $pCompoundDB, $PrintDetail);

	$Results[$c]->{MinE}            = $MinE;
	$Results[$c]->{pMinComposition} = $pMinComposition;
	$c++;
}

&mprint($check, "\n");
&mprint($check, "Minimum energy for [$Reaction]\n");
my $NormalizedComposition;
for(my $i = 0 ; $i < $c ; $i++) {
	my $r = $Results[$i];
	my $MinE            = $r->{MinE};
	my $pMinComposition = $r->{pMinComposition};
	my $pnFinalElement  = $R->nElementByArray(\@PossibleCompounds, $pMinComposition);
	my $FinalProducts   = $R->BuildReactionFromArrays(\@PossibleCompounds, $pMinComposition, undef, undef, 1);
	$NormalizedComposition = $R->NormalizeComopositionByArray(\@PossibleCompounds, $pMinComposition, "%4.2f");
	&mprint($check, "  Iter$i:\n");
	&mprintf($check, "    Stable composition: %10.4f eV: %s\n", $MinE, $FinalProducts);
	my $Eform = $MinE - $EReaction;
	&mprintf($check, "                        Eform = %10.4f eV\n", $Eform);
#	&mprintf($check, "    Last composition  : %10.4f eV: %s\n", $EReaction, $NormalizedComposition);
	$log->printf("  Iter$i: %10.4f eV: %s\n", $MinE, $FinalProducts);
	$log->print("    Last compositiion: ", $R->NormalizeComopositionByArray(\@PossibleCompounds, $pMinComposition, "%4.2f"), "\n");
}

&mprint($check,  "  Last:\n");
&mprintf($check, "    Last composition  : %10.4f eV: %s\n", $EReaction, $Reaction);
&mprintf($check, "                        %s\n", $NormalizedComposition);
#&mprint($check, "      Etot for [$Reaction]: $EReaction eV\n");

$log->print("\n");
$log->Close();
$check->Close();
exit;

#==============================================
# 最小化関数のサブルーチン
# MinimizationFunc    : 
# PDLDisplaySimplex: PDL::Opt::Simplex用のDisplay関数
#==============================================
sub MinimizationFunc {
	my ($pComposition, $iPrintLevel) = @_;
	$iPrintLevel = 1 if($iPrintLevel > 1);

	$iIteration++;

	$R->RepairComposition($pComposition, $pReagentElements, \@PossibleCompounds, $pCompoundDB, $PrintDetail);
	my $E = $R->CalTotalEnergyByArray(\@PossibleCompounds, $pComposition, $PrintDetail);

	my $Products = $R->BuildReactionFromArrays(\@PossibleCompounds, $pComposition, undef, undef, 1);
my $s = $R->NormalizeComopositionByArray(\@PossibleCompounds, $pComposition, "%4.2f");
#printf("E=%8.4f eV: %s [%s]\n", $E, $Products, $s);
	if($iIteration % $iPrintInterval == 0) {
		printf("E=%8.4f eV: %s [%s]\n", $E, $Products, $s) if($iPrintLevel);
	}
	return $E;
}

sub PDLDisplayFunc {
	my ($piddle) = @_;
#print "Display: piddle=$piddle\n";
}

sub mprint
{
	my ($out, @a) = @_;
	print(@a);
	$out->print(@a);
}

sub mprintf
{
	my ($out, $format, @a) = @_;
	printf($format, @a);
	$out->printf($format, @a);
}