#!/usr/bin/perl

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

use strict;

use Utils;
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 = "TotalEnergy-cubic-a.csv";
$SourceCSV    = $ARGV[0] if(defined $ARGV[0]);

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

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 ($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 $OutCSV    = Deps::MakePath("$drive$directory", "${filebody}-Out.csv", 0);
my $FitCSV    = Deps::MakePath("$drive$directory", "${filebody}-Fit.csv", 0);

my $pInf  = &ReadCSV($SourceCSV);
my $pPcal = &CalPressure($SourceCSV, $pInf);
my $out = JFile->new($OutCSV, "w") or "$!: Can not write to [$OutCSV].\n";
$out->print("Pcal(GPa),P(GPa),a(A),b(A),c(A),V(A^3),Etot(eV),H/Z(eV),Z\n");
for(my $i = 0 ; $i < @$pInf ; $i++) {
	my $V = $pInf->[$i]{V};
	next if($V < $FitVMin or $FitVMax < $V);

	my $E = $pInf->[$i]{E};
	my $P = (defined $pInf->[$i]{P})? $pInf->[$i]{P} : $pPcal->[$i];
	my $H = ($E + $P * 1.0e9 * $V * 1.0e-30 / $e) / $pInf->[$i]{Z};
	$out->print("$pPcal->[$i],$pInf->[$i]{P},$pInf->[$i]{a},$pInf->[$i]{b},$pInf->[$i]{c},$pInf->[$i]{V},$E,$H,$pInf->[$i]{Z}\n");
}
$out->Close();

my ($Emin, $Vmin, $B, $B1, $MinVal) = &CalEOS($SourceCSV, $pInf);

exit;

#========================================================================
# Subroutines
#========================================================================
sub ReadCSV
{
	my ($SourceCSV) = @_;

	my $csv = new CSV;
	$csv->Read($SourceCSV) or die "$!: Can not read [$SourceCSV].\n";
	$csv->Close();
	my $pZ = $csv->GetXData("Z");
	my $pa = $csv->GetXData("a");
	my $pb = $csv->GetXData("b");
	my $pc = $csv->GetXData("c");
	my $pV = $csv->GetXData("V.*");
	my $pP = $csv->GetXData("PStress.*");
	my $pE = $csv->GetXData("Etot.*");
	my $nData = @$pV;

print "Source: [$SourceCSV]\n";
print "Out   : [$OutCSV]\n";
print "pa=$pa\n";
print "pb=$pb\n";
print "pc=$pc\n";
print "pV=$pV\n";
print "pP=$pP\n";
print "pE=$pE\n";
print "nData=$nData\n";

	my @Inf;
	for(my $i = 0 ; $i < $nData ; $i++) {
		$Inf[$i] = {
			Z => $pZ->[$i],
			a => $pa->[$i],
			b => $pb->[$i],
			c => $pc->[$i],
			V => $pV->[$i],
			P => $pP->[$i] * 0.1, # kbar => GPa
			E => $pE->[$i],
			};
	}
	my @sorted = sort { return $a->{P} <=> $b->{P} if(defined $a->{P});
						return $b->{V} <=> $a->{V} } @Inf;
#for(my $i = 0 ; $i < @sorted ; $i++) {
#print "i=$i: P=$sorted[$i]{P}, V=$sorted[$i]{V}\n";
#}
#exit;

	return \@sorted;
}

sub CalPressure
{
	my ($SourceCSV, $pInf) = @_;

	my $nData = @$pInf;
	my $nl   = $nData - 1;
	my ($pE, $pV);
	for(my $i = 0 ; $i < $nData ; $i++) {
		$pE->[$i] = $pInf->[$i]{E};
		$pV->[$i] = $pInf->[$i]{V};
	}

	my @Pcal;
	$Pcal[0]   = ($pE->[1]   - $pE->[0])     / ($pV->[1]   - $pV->[0]);
	$Pcal[$nl] = ($pE->[$nl] - $pE->[$nl-1]) / ($pV->[$nl] - $pV->[$nl-1]);
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		$Pcal[$i] = 0.5 * (
						  ($pE->[$i+1] - $pE->[$i]) / ($pV->[$i+1] - $pV->[$i])
						+ ($pE->[$i] - $pE->[$i-1]) / ($pV->[$i] - $pV->[$i-1])
						);
	}
	for(my $i = 0 ; $i < $nData ; $i++) {
		$Pcal[$i] *= -$e * 1.0e30 * 1.0e-9; #GPa
	}

	return \@Pcal;
}

sub CalEOS
{
	my ($SourceCSV, $pInf) = @_;

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

	my $nData = @$pInf;
	my $nl   = $nData - 1;
	my ($pE, $pV);
	for(my $i = 0 ; $i < $nData ; $i++) {
		$pE->[$i] = $pInf->[$i]{E};
		$pV->[$i] = $pInf->[$i]{V};
	}

	my $txt = JFile->new($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($pInf, @_); },
					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($pInf, @_); },
					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";
	
	return ($Emin, $Vmin, $B, $B1, $MinVal);}

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

	my $nData = @$pInf;
	my $nl   = $nData - 1;
	my ($pE, $pV);
	for(my $i = 0 ; $i < $nData ; $i++) {
		$pE->[$i] = $pInf->[$i]{E};
		$pV->[$i] = $pInf->[$i]{V};
	}

	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 ($pInf, $pVars, $iPrintLevel) = @_;
	
	my $nData = @$pInf;
	my $nl   = $nData - 1;
	my ($pE, $pV);
	for(my $i = 0 ; $i < $nData ; $i++) {
		$pE->[$i] = $pInf->[$i]{E};
		$pV->[$i] = $pInf->[$i]{V};
	}

	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;
}