#!/usr/bin/perl

use strict;

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

use IniFile;
use CSV;
use Sci::Optimize;

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

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

my $ObsCSVPath  = "20110217-5nmonSi.csv";
#my $ObsCSVPath  = "20110217-10nmonSi.csv";
my $SiCSVPath   = "20110217onlySi.csv";
my $IGZOCSVPath = "20110217-60nm.csv";
my $OutputCSV   = "Output.csv";

print "Read [$ObsCSVPath].\n";
my $ObsCSV  = new CSV();
$ObsCSV->Read($ObsCSVPath, 1)   or die "$!: Can not read [$ObsCSVPath].\n";g:
&ReverseX($ObsCSV);
my $pXObs = $ObsCSV->GetXData("E.*", ".*eV.*");
my $pYObs = $ObsCSV->GetYData(".*signal.*", "I.*");
$ObsCSV->CalMinMax();
print "pXObs=$pXObs  pY=$pYObs\n";

print "Read [$SiCSVPath].\n";
my $SiCSV   = new CSV();
$SiCSV->Read($SiCSVPath, 1)     or die "$!: Can not read [$SiCSVPath].\n";
&ReverseX($SiCSV);
my $pXSi = $SiCSV->GetXData("E.*", ".*eV.*");
my $pYSi = $SiCSV->GetYData(".*signal.*", "I.*");
print "pXSi=$pXSi  pY=$pYSi\n";
$SiCSV->CalMinMax();

print "Read [$IGZOCSVPath].\n";
my $IGZOCSV = new CSV();
$IGZOCSV->Read($IGZOCSVPath, 1) or die "$!: Can not read [$IGZOCSVPath].\n";
&ReverseX($IGZOCSV);
my $pXIGZO = $IGZOCSV->GetXData("E.*", ".*eV.*");
my $pYIGZO = $IGZOCSV->GetYData(".*signal.*", "I.*");
print "pXIGZO=$pXIGZO  pY=$pYIGZO\n";
$IGZOCSV->CalMinMax();

my $nData         = $ObsCSV->nData();
my ($xmin, $xmax) = $ObsCSV->GetXMinMax();
my $xstep = ($xmax - $xmin) / ($nData - 1);
print "$nData: ($xmin - $xmax, $xstep)\n";

my $optimize = new Optimize;
print "use $Method\n";

#=========================================================
# 最適化変数設定
#=========================================================
	my (@X, @Y);
	my $c = 0;
	for(my $i = 0 ; $i < @$pXObs ; $i++) {
		my $x = $pXObs->[$i];
#		next if($x < $XMin or $XMax < $x);
#		next if($nSkipData > 0 and $i % $nSkipData != 0);
		$X[$c] = $x;
		$Y[$c] = $pYObs->[$i];
		$c++;
	}
	my (@pSpectrumData) = ($SiCSV, $IGZOCSV);

	my $dini   = new IniFile(undef, 0, 0);
	my $nSpectra = @pSpectrumData;
	&AddParameters($dini, "BG",   0.0, 1, 1000.0, 0.0, '');
	for(my $i = 0 ; $i < $nSpectra ; $i++) {
		&AddParameters($dini, "k$i",  1.0, 1, 0.1, 0.0, '');
		&AddParameters($dini, "dE$i", 0.0, 1, 0.1,  '', '');
	}

#=========================================================
# 最適化オブジェクト作成
#=========================================================
	$optimize->AddParameters(
			"BG",  \$dini->{BG},
			$dini->{BGcheck}, $dini->{BGscale}, $dini->{BGmin}, $dini->{BGmax},
			sub {
				$dini->{BG} = $_[0];
			},
		);
	for(my $i = 0 ; $i < $nSpectra ; $i++) {
		my $ii = $i;
		$optimize->AddParameters(
			"k$ii",  \$dini->{"k$ii"},
			$dini->{"k${ii}check"}, $dini->{"k${ii}scale"}, $dini->{"k${ii}min"}, $dini->{"k${ii}max"},
			sub {
				$dini->{"k$ii"} = $_[0];
			},

			"dE$ii",  \$dini->{"dE$ii"},
			$dini->{"dE${ii}check"}, $dini->{"dE${ii}scale"}, $dini->{"dE${ii}min"}, $dini->{"dE${ii}max"},
			sub {
				$dini->{"dE$ii"} = $_[0];
			},
		);
	}

#=========================================================
# 最適化の実行
#=========================================================
	$optimize->SetnS2Calculation(0);
	my ($OptVars, $MinVal) = $optimize->Optimize(
				$Method, undef, undef, undef,
				$EPS, $nMaxIter, $iPrintLevel,
				sub { &CalS2(\@X, \@Y, \@pSpectrumData, @_); },
				undef,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
	print "\nOptimized at:\n";

	$optimize->RecoverParameters($OptVars);
	$optimize->PrintParameters(1, $OptVars, $MinVal, 1);
exit;

#==============================================
# サブルーチン
#==============================================
sub AddParameters
{
	my ($dini, @args) = @_;
	for(my $i = 0 ; $i < @args ; $i += 6) {
		my $key = $args[$i];
		$dini->{$key}          = $args[$i+1] if(!defined $dini->{$key});
		$dini->{"${key}check"} = $args[$i+2] if(!defined $dini->{"${key}check"});
		$dini->{"${key}scale"} = $args[$i+3] if(!defined $dini->{"${key}scale"});
		$dini->{"${key}min"}   = $args[$i+4] if(!defined $dini->{"${key}min"});
		$dini->{"${key}max"}   = $args[$i+5] if(!defined $dini->{"${key}max"});
#print "AddParameters: $key: $dini->{$key}\n";
	}
}

sub ReverseX
{
	my ($csv) = @_;
return $csv->ReverseData();
	
	my $n          = $csv->nData();
	my $pDataArray = $csv->DataArray();
	for(my $i = 0 ; $i < @$pDataArray ; $i++) {
		my $p = $pDataArray->[$i];
		for(my $j = 0 ; $j < $n / 2 ; $j++) {
			($p->[$j], $p->[$n-1-$j]) = ($p->[$n-1-$j], $p->[$j]);
		}
	}
}

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

	my $BG = $ini->{BG};
	my (@k, @dE);
	my @pCalY;
	my @YTotal;
	for(my $j = 0 ; $j < @$pX ; $j++) {
		$YTotal[$j] = 0.0;
	}

	for(my $i = 0 ; $i < $nSpectra ; $i++) {
		my $k     = $ini->{"k$i"};
		my $dE    = $ini->{"dE$i"};
		my $pData = $ppSpectrumData->[$i];
		next if(!defined $pData);

		my @YCal;
		for(my $j = 0 ; $j < @$pX ; $j++) {
			my $x        = $pX->[$j];
			my $y        = $k * $pData->YVal($x + $dE);
			$YCal[$j]    = $y;
			$YTotal[$j] += $y;
		}
		$pCalY[$i] = \@YCal;
	}
	return (\@YTotal, \@pCalY);
}

sub CalS2
{
	my ($pX, $pY, $ppSpectrumData, $pVars, $iPrintLevel) = @_;
#print "pVars=", join(',', @$pVars), "\n";

	my $ini = $dini;
	$optimize->RecoverParameters($pVars);

	my %p;
	my $idx = 0;
	$p{BG} = $pVars->[$idx++];
	for(my $i = 0 ; $i < $nSpectra ; $i++) {
		$p{"k$i"}  = $pVars->[$idx++];
		$p{"dE$i"} = $pVars->[$idx++];
	}
	my ($pYTotal, $ppCalY) = &CalTotalSpectrum($pX, $ppSpectrumData, \%p);

	my $S2 = 0.0;
	my $c = 0;
	for(my $i = 0 ; $i < @$pX ; $i++) {
		my $x = $pX->[$i];
		my $d = $pYTotal->[$i] - $pY->[$i];
		$S2 += $d * $d;
		$c++;
	}
	$S2 = sqrt($S2 / $c) if($c > 0);

	$optimize->IncrementnS2();
	my $nS2         = $optimize->nS2Calculation();
	print "S2[$nS2] = $S2 / ";
	$S2 += $optimize->CalPenalty(1.0e10, $pVars, 1);
	print "$S2\n";
	$optimize->PrintParameters(1, $pVars, $S2, 0);

	return $S2;
}