#!/usr/bin/perl

use strict;

use lib 'd:/Programs/Perl/lib';
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 = 100;
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 @Guess = (0.5,  0.0, 0.5,  0.0);
my @Scale = (0.1, 0.05, 0.1, 0.05);
my @OptId = (  1,    1,   1,    1);
my $Opt = new Optimize;
print "use $Method\n";

my ($OptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
				$EPS, $nMaxIter, $iPrintLevel,
				\&MinimizationFunc, \&PDLDisplayFunc,
				sub { Optimize::BuildDifferentialMatrixes(@_); },
				);
print "Optimized at (",join(',',@{$OptVars}),")=$MinVal\n";

	my ($kSi, $dESi, $kIGZO, $dEIGZO) = @$OptVars;
	my $out = new JFile;
	$out->Open($OutputCSV, "w") or die "$!: Can not write to [$OutputCSV].\n";
	$out->print("E(eV),Obs,Calc,Si,IGZO\n");
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $E     = $xmin + $i * $xstep;
		my $y     = $ObsCSV->YVal( $pXObs,  $pYObs,  $E);
		my $ySi   = $SiCSV->YVal(  $pXSi,   $pYSi  , $E-$dESi);
		my $yIGZO = $IGZOCSV->YVal($pXIGZO, $pYIGZO, $E-$dEIGZO);
		my $sum   = $kSi * $ySi + $kIGZO * $yIGZO;
		$out->print("$E,$y,$sum,$ySi,$yIGZO\n");
	}
	$out->Close();
	
	print "Save fitting results to [$OutputCSV].\n";

exit;

#==============================================
# サブルーチン
#==============================================
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 MinimizationFunc {
	my ($pVars, $iPrintLevel) = @_;
	my ($kSi, $dESi, $kIGZO, $dEIGZO) = @$pVars;

	my $S = 0.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $E     = $xmin + $i * $xstep;
		my $y     = $ObsCSV->YVal( $pXObs,  $pYObs,  $E);
		my $ySi   = $SiCSV->YVal(  $pXSi,   $pYSi,   $E+$dESi);
		my $yIGZO = $IGZOCSV->YVal($pXIGZO, $pYIGZO, $E+$dEIGZO);
#my $min = $pXObs->[0];
#my $max = $pXObs->[$nData-1];
#print "$i: $E  $y ($min - $max)\n";

		my $d = $y - $kSi * $ySi - $kIGZO * $yIGZO;
		$S += $d*$d;
	}
	$S = sqrt($S/($nData-1));
printf("kSi,dESi,kIGZO,dEIGZO=%6.3f\t%6.3f\t%6.3f\t%6.3f\t(%g)\n", 
		$kSi, $dESi, $kIGZO, $dEIGZO, $S) if($iPrintLevel >= 2);
	return $S;
}

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