#!/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 Sci::ChemicalReaction;
use Sci::Optimize;

my $Opt = new Optimize;

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

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

print "use $Method\n";

my $PrintDetail = 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 '')? "La+Fe+O+As+0.1As+0.1O+0.1H" : $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 @Products = (
		"2SrCuFSe",
		"MnP",
		"BaO2H2",
		"BaO+H2O",
		"EuCuFSe+0.1Cu",
		"1/18Eu18Cu17F18Se18+0.1555556Cu",
		);
my $BaseDBDir = "Databases";
my @DBFiles   = ("AeCuFCh-High.csv",
			"AeCuFCh.csv", "LaMnOPn.csv", "BaFe2As2.csv", "SrFe2As2.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",
			);

#=============================================
# Create ChemicalReaction object
#=============================================
my $R = new ChemicalReaction(\@DBFiles, $BaseDBDir);
if(!$R->DB()) {
	print "Error: $!: Can not read a DB.\n";
	exit;
}
my $pCompoundDB   = $R->CompoundDB();

if(0) {
my $s = "Sr[(OH)2(OCl4)3]2Br";
my ($pe, $pn) = $R->CompoundToElements($s, 1);
print "s: $s\n";
for(my $i = 0 ; $i < @$pe ; $i++) {
	print "$pe->[$i]: $pn->[$i]\n";
}
exit;
}

#=============================================
# 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);
	print "$e: n=$n: Z=$pAtomDB->{Z}: M=$pAtomDB->{mass}\n";
}

my @PossibleCompounds = $R->ExtractPossibleCompounds();#$pElements, [keys %Compound]);
print "To be considered: Atom\n";
foreach my $c (sort @PossibleCompounds) {
	next if($pCompoundDB->{$c}->{Type} ne 'Atom');
	printf "  %-2s %6.3f eV", $c, $pCompoundDB->{$c}->{Emol};
	print  " ($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');
	printf "  %-2s %6.3f eV", $c, $pCompoundDB->{$c}->{Emol};
	print  " $pCompoundDB->{$c}->{Type} ($pCompoundDB->{$c}->{Compound} in $pCompoundDB->{$c}->{Ref})\n";
}
print "To be considered: Compound\n";
foreach my $c (sort @PossibleCompounds) {
	next if($pCompoundDB->{$c}->{Type} ne 'Compound');
	printf "  %s %10.3f eV", $c, $pCompoundDB->{$c}->{Emol};
	print  " $pCompoundDB->{$c}->{Type} ($pCompoundDB->{$c}->{Compound} in $pCompoundDB->{$c}->{Ref})\n";
}

my $EReaction = $R->CalTotalEnergy($Reaction);
print "Reaction: $Reaction\n";
print "  E(Reaction)=$EReaction eV\n";

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

my @Results;
my $c = 0;
for(my $l = -2 ; $l <= -2 ; $l++) {
#for(my $l = -2 ; $l < @PossibleCompounds ; $l++) {
	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, $PrintDetail);
		$k = $KScale0;
	}
	elsif($l == -1) {
		for(my $i = 0 ; $i < @Guess ; $i++) {
			$pComposition->[$i] = $MaxNInReaction;
		}
		$pComposition = $R->RepairComposition($pComposition, $pReagentElements, \@PossibleCompounds, $PrintDetail);
		$k = $KScale0;
	}
	else {
		$pComposition = $R->MakePossibleComposition($pReagentElements, \@PossibleCompounds, $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, $PrintDetail);

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

print "\n";
print "Minimum energy for [$Reaction]\n";
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);
	printf "  %10.4f eV: %s\n", $MinE, $FinalProducts;
	print "  ";
	foreach my $e (sort keys %$pnFinalElement) {
		print "$e$pnFinalElement->{$e}";
	}
	print "\n";
}
print "Etot for [$Reaction]: $EReaction eV\n";
for(my $i = 0 ; $i < @Products ; $i++) {
	my $E = $R->CalTotalEnergy($Products[$i], 0);
	print "Etot for [$Products[$i]]: $E eV\n";
}

exit;

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

	$iIteration++;

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

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

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