#!/usr/bin/perl

use lib 'c:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';

use strict;

use Sci::Algorism;
use JFile;
use CSV;

use Sci qw($pi $kB $R $NA $c $e $e0 $u0 $me $mp $mn $h $hbar $torad $todeg);
use Sci::Optimize;

my $SourceCSV = "NaCl-NaCl.csv";
$SourceCSV    = $ARGV[0] if(defined $ARGV[0]);

#my $EOSModel = "Murnaghan";
my $EOSModel = "Birch-Murnaghan";

my ($FitVMin, $FitVMax) = (0, 1200);
my $nSimplexIter = 1;

my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($SourceCSV);
my $OutTXT    = Deps::MakePath("$drive$directory", "${filebody}-Out.txt", 0);
my $FitCSV    = Deps::MakePath("$drive$directory", "${filebody}-Fit.csv", 0);

my $Method = "Amoeba::Simplex";
#my $Method = "PDL::Simplex";
#my $Method = "ModifiedNewton";
#my $Method = "LinearOptimization"; # 線形パラメータのみ。
my $EPS         = 1.0e-7;
my $nMaxIter    = 3000;
my $iPrintLevel = 2;

my $csv = new CSV;
$csv->Read($SourceCSV) or die "$!: Can not read [$SourceCSV].\n";
my $pV = $csv->GetXData("V.*");
print "pV=$pV\n";
my $pE = $csv->GetXData("Etot.*");
print "pE=$pE\n";
my $nData = @$pV;
print "nData=$nData\n";

my $txt = new JFile;
$txt->Open($OutTXT, "w") or die "$!: Can not write to [$OutTXT].\n";
$txt->mprint("Model    : $EOSModel\n");
$txt->mprint("Source   : $SourceCSV\n");
$txt->mprint("FitCSV   : $FitCSV\n");

my ($Vmin, $Emin) = (0, 1.0e10);
for(my $i = 0 ; $i < $nData ; $i++) {
	if($Emin > $pE->[$i]) {
		$Vmin = $pV->[$i];
		$Emin = $pE->[$i];
	}
}
$txt->mprint("Vmin = $Vmin @ E=$Emin\n");

my $m = 5;
my $Opt = new Optimize;
my ($OptVars, $MinVal) = $Opt->Optimize("mlsq", $pV, $pE, $m, $iPrintLevel);

my @Guess     = @$OptVars;
my @OptId = (   1,   1,   1,   1);
my @Scale = (   10,   10,   10,   10);
for(my $i = 0 ; $i < $nSimplexIter ; $i++) {
	my $Opt = new Optimize;
	($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				sub { &MinimizationMLSQFunc(@_); },
				sub { },
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
	@Guess = @$OptVars;
	print "Optimized at (",join(',',@{$OptVars}),")=$MinVal\n";
}
my @c = @$OptVars;
$txt->mprint("\n");
$txt->mprint("Optimized parameters\n");
for(my $i = 0 ; $i <= $m ; $i++) {
	$txt->mprint("c$i: $c[$i]\n");
}
$txt->mprint("MinS: $MinVal\n");
$txt->mprint("\n");

my $B = 2.0 * $c[2] + 6.0 * $c[3] * $Emin + 12.0 * $c[4] * $Emin * $Emin;
$txt->mprint("B=$B\n");

my $B1 = 3.5;
my @Guess = ($Emin, $Vmin, $B, $B1);
my @OptId = (    1,     1,  1,   1);
my @Scale = (   10,    10,  1,   1);
for(my $i = 0 ; $i < $nSimplexIter ; $i++) {
	my $Opt = new Optimize;
	($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				sub { &MinimizationEOSFunc(@_); },
				sub { },
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
	@Guess = @$OptVars;
	print "Optimized at (",join(',',@{$OptVars}),")=$MinVal\n";
}
($Emin, $Vmin, $B, $B1) = @Guess;

$txt->mprint("\n");
$txt->mprint("Optimized parameters\n");
$txt->mprint("E0 : $Emin eV\n");
$txt->mprint("V0 : $Vmin A^3\n");
$txt->mprint("B0 : $B, ", $B * $e * 1.0e30 * 1.0e-9, " GPa\n");
$txt->mprint("B0': $B1\n");
$txt->mprint("Min: $MinVal\n");
$txt->mprint("\n");

my $fit = new JFile;
$fit->Open($FitCSV, "w") or "$!: Can not write to [$FitCSV].\n";
$fit->print("P(GPa),V(A^3),Emlsq(eV),EEOS(eV),EDFT(eV)\n");
for(my $i = 0 ; $i < $nData ; $i++) {
	my $V = $pV->[$i];
	next if($V < $FitVMin or $FitVMax < $V);

	my $E = $pE->[$i];
	my $EcalMLSQ = &CalEMLSQ($V, @c);
	my $EcalEOS  = &CalEEOS($V, $Emin, $Vmin, $B, $B1);
	# CalPEOS: P = dE / dV [eV / A3]
	my $PEOS     = &CalPEOS($V, $Emin, $Vmin, $B, $B1) * $e * 1.0e30 * 1.0e-9; # GPa
	print("$V\t$PEOS\t$EcalMLSQ\t$EcalEOS\t$E\n");
	$fit->print("$PEOS,$V,$EcalMLSQ,$EcalEOS,$E\n");
}
$fit->Close();
$txt->Close();

print "Vmin = $Vmin @ E=$Emin\n";

exit;

#==============================================
# 最小化関数のサブルーチン
# MinimizationFunc    : 
# PDLDisplaySimplex: PDL::Opt::Simplex用のDisplay関数
#==============================================
sub MinimizationEOSFunc {
	my ($pVars, $iPrintLevel) = @_;
	
	my ($E0, $V0, $B0, $B1) = @$pVars;
	my $S = 0.0;
	my $c = 0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $V = $pV->[$i];
		next if($V < $FitVMin or $FitVMax < $V);

		my $E = $pE->[$i];
		my $Ecal = &CalEEOS($V, $E0, $V0, $B0, $B1);
		$S += ($Ecal - $E) * ($Ecal - $E);
	
		$c++;
	}
	$S += 1.0e10 if($B1 <= 1.0);

if($iPrintLevel >= 2) {
	printf("E0,V0,B0,B1=%5.3g, %5.3f, %5.3g, %5.3g ", $E0, $V0, $B0, $B1, $S);
	printf("(%g)\n", $S);
}

	return $S / $c;
}

sub CalEEOS
{
	my ($V, $E0, $V0, $B0, $B1) = @_;

	if($EOSModel eq "Murnaghan") {
		return $E0 + $B0 * $V / $B1 * ( ($V0/$V)**$B1 / ($B1 - 1.0) + 1.0) - $B0 * $V0 / ($B1 - 1.0);
#		return $E0 + $B0 * $V / $B1 * ( ($V0/$V)**$B1 / ($B1 - 1.0) + 1.0) + $B0 * $B0 / ($B1 - 1.0);
	}
	else {
		my $rV   = $V0 / $V;
		my $rV23 = $rV**(2.0/3.0);
		return $E0 + 9.0/16.0*$V0*$B0 * (
					($rV23 - 1.0)**3.0 * $B1 + ($rV23 - 1.0)**2.0 * (6.0 - 4.0*$rV23)
					);
	}
}

sub CalPEOS
{
	my ($V, $E0, $V0, $B0, $B1) = @_;

	if($EOSModel eq "Murnaghan") {
		return $B0 / $B1 * ( ($V0/$V)**$B1 - 1.0);
	}
	else {
		my $rV   = $V0 / $V;
		my $rV23 = $rV**(2.0/3.0);
		my $rV53 = $rV**(5.0/3.0);
		my $rV73 = $rV**(7.0/3.0);
		return 3.0/2.0*$B0 * ( $rV73 - $rV53) * (1.0 + 3.0/4.0*($B1 - 4.0) * ($rV23 - 1.0) );
	}
}


sub MinimizationMLSQFunc {
	my ($pVars, $iPrintLevel) = @_;
	
	my @c = @$pVars;
	my $S = 0.0;
	my $c = 0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $V = $pV->[$i];
		next if($V < $FitVMin or $FitVMax < $V);

		my $E = $pE->[$i];
		my $Ecal = &CalEMLSQ($V, @c);
		$S += ($Ecal - $E) * ($Ecal - $E);
	
		$c++;
	}
if($iPrintLevel >= 2) {
	printf("c=%5.3g, %5.3f, %5.3g, %5.3g ", @c, $S);
	printf("(%g)\n", $S);
}

	return $S / $c;
}

sub CalEMLSQ
{
	my ($V, @c) = @_;
	my $ret = $c[0];
	my $Vn = $V;
	for(my $i = 1 ; $i < @c ; $i++) {
		$ret += $c[$i] * $Vn;
		$Vn *= $V;
	}
	return $ret;
}