#!/usr/bin/perl

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

use strict;

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

my $SourceCSV    = "NaCl-NaCl.csv";
#my $SourceCSV    = "NaCl-CsCl.csv";
#my $SourceCSV    = "CaO-CsCl.csv";
#my $SourceCSV    = "CaO-NaCl.csv";
my ($FitVMin, $FitVMax) = (0, 70);
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("Emol.*");
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("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\n");
$txt->mprint("V0 : $Vmin\n");
$txt->mprint("B0 : $B\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("V(A^3),P(),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);
	my $PEOS     = &CalPEOS($V, $Emin, $Vmin, $B, $B1);
	print("$V\t$PEOS\t$EcalMLSQ\t$EcalEOS\t$E\n");
	$fit->print("$V,$PEOS,$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) = @_;
	return $E0 + $B0 * $V / $B1 * ( ($V0/$V)**$B1 / ($B1 - 1.0) + 1.0) + $B0 * $B0 / ($B1 - 1.0);
}

sub CalPEOS
{
	my ($V, $E0, $V0, $B0, $B1) = @_;
	return $B0 / $B1 * ( ($V0/$V)**$B1 - 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;
}