#!/usr/bin/perl

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

use strict;
use Sci::Optimize;

my $Opt = new Optimize;

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

my $EPS = 1.0e-7;
my $nMaxIter = 100;
my @Guess = (1, 1, 1, 1); #線形最小自乗法の場合、初期値は意味がない。@Guessの要素数が変数の数になる
my @OptId = (1, 1, 1, 1);
my @Scale = (1, 1, 1, 1);
my $iPrintLevel = 0;

print "use $Method\n";

my $nData = 30;

#線形最小自乗法の場合、\@Scale, $EPS, $nMaxIter, を渡す必要はない
#必要なければ、PDLDisplayFunc以降も指定する必要はない
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 $OutFile = "a.csv";
open(OUT, ">$OutFile") or die "Can not write to [$OutFile].\n";
print OUT "x,yobs,ycal\n";
for(my $i = 0 ; $i < $nData ; $i++) {
	my $x = IdxToX($i);
	my $yobs = YObs($x);
	my $ycal = YCal($x, @$OptVars);
	print OUT "$x,$yobs,$ycal\n";
}
close(OUT);
exit;

#==============================================
# 最小化関数のサブルーチン
# MinimizationFunc    : 
# PDLDisplaySimplex: PDL::Opt::Simplex用のDisplay関数
#==============================================
sub IdxToX
{
	my ($i) = @_;
	return $i * 0.1;
}

sub YObs
{
	my ($x) = @_;
	my $yobs = exp(-$x*$x);
	return $yobs;
}

sub YCal
{
	my ($x, $c0, $c1, $c2, $c3) = @_;
	my $ycal = $c0 + $c1 * $x + $c2 * $x*$x + $c3 * $x*$x*$x;
	return $ycal;
}

sub MinimizationFunc {
	my ($pVars, $iPrintLevel) = @_;
	my ($c0, $c1, $c2, $c3) = @$pVars;

	my $S = 0.0;
	for(my $i = 0 ; $i < $nData ; $i++) {
		my $x = IdxToX($i);
		my $yobs = YObs($x);
		my $ycal = YCal($x, @$pVars);
		my $d = $yobs - $ycal;
		$S +=  $d*$d;
	}
#printf("c=%6.3f\t%6.3f\t%6.3f\t%6.3f\t(%g)\n", $c0, $c1, $c2, $c3, $S) if($iPrintLevel >= 2);
	return $S / $nData;
}

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