#!/usr/bin/perl

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

use strict;

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

my $SourceCSV = "d-dependence.csv";
my $FitCSV    = "d-dependence-Fit.csv";

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

my $EPS = 1.0e-7;
my $nMaxIter = 1000;
my $iPrintLevel = 2;

my $csv = new CSV;
if($csv->Read($SourceCSV)) {
	print("  Read [$SourceCSV]\n");
}
else {
	print("  Error!!: Can not read [$SourceCSV]\n");
	exit;
}
my $pLabelArray = $csv->LabelArray();
print "Label: ", join(',', @$pLabelArray), "\n";
my $pDataArray  = $csv->DataArray();
my $pd = $csv->GetXData("d.*");
#my $pI = $csv->GetYData("ISi-O.*");
my $pI = $csv->GetYData("ISi1s.*");
my $nData = @$pd;
print "nData = $nData\n";

my $delta = 7.0; # nm
my $I0    = $pI->[0] / $delta;

my @Guess = ( $delta,        $I0);
my @OptId = (      1,          1);
my @Scale = ($I0*1, $delta*1);
print "use $Method\n";

my $Opt = new Optimize;

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

my $out = new JFile;
$out->Open($FitCSV, "w") or die "$!: Can not write to [$FitCSV]\n";
$out->print("d(nm),Iobs,Ical\n");
for(my $i = 0 ; $i < $nData ; $i++) {
	my $d = $pd->[$i];
	my $Ical = &Cal($d, $delta, $I0);
	$out->print("$d,$pI->[$i],$Ical\n");
}
$out->print("delta=,$delta,nm\n");
$out->print("I0=,$I0\n");
$out->Close();

exit;

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

	my $S = 0.0;
	my ($delta, $I0) = @$pVars;
	for(my $i = 0 ; $i < $nData ; $i++) {
		next if(!defined $pI->[$i] or $pI->[$i] eq '');
		my $d = $pd->[$i];
		my $Ical = &Cal($d, $delta, $I0);
		$S += ($Ical - $pI->[$i])**2;
	}
	$S += 1.0e5 if($I0 < 0.0);
	$S += 1.0e5 if($delta < 0.0);
printf("delta,I0==%6.3f\t%6.3f\t(%g)\n", $delta, $I0, $S) if($iPrintLevel >= 2);
	return $S;
}

sub Cal
{
	my ($d, $delta, $I0) = @_;
	return $I0 * $delta * exp(-$d / $delta);
}

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