end0tknr's kipple - 新web写経開発

http://d.hatena.ne.jp/end0tknr/ から移転しました

Math::GSL ( GNU Scientific Library ) for perl による線形補間とスプライン補間

linuxで科学技術計算と言えば、 GSL( GNU Scientific Library )を思い浮かべますが、GSLのwrapper である Math::GSL for perl を見かけたので、リファレンス・マニュアル(日本語版のpdfもあります)にあるサンプルコードを参考に、線形補間とスプライン補間を写経。
Math::GSL::Interp - search.cpan.org
GSL - GNU Scientific Library - GNU Project - Free Software Foundation (FSF)
GSL reference manual Japanese translation

#!/usr/local/bin/perl
use strict;
use warnings;
use utf8;
use Encode;
use Data::Dumper;

# use Math::GSL::Heapsort qw /:all/;
use Math::GSL::Interp qw/:all/;
# use Math::GSL::Sort qw/:all/;
# use Math::GSL::Statistics qw /:all/;

# http://search.cpan.org/~ambs/Math-GSL/lib/Math/GSL/Interp.pm
# http://www.gnu.org/software/gsl/
# http://cbrc3.cbrc.jp/~tominaga/translations/gsl/index.html (pdf manual)

main();


sub main {


    #sin(), cos()により、オジリナルの擬似曲線を用意
    my @x;
    my @y;
    print "GENERATED ORIGINAL CURVE by sin() & cos()\n";
    for (my $i=0; $i<10; $i++) {
        $x[$i] = sprintf("%.2f",  $i + 0.5 * sin($i));
        $y[$i] = sprintf("%.2f",  $i + cos($i * $i) );
        print "X=$x[$i] , Y=$y[$i]\n";
    }


    my $interp_vals = {};
    my @interp_x_vals = qw/0 1 1.5  3.5  5  7.2  8.6  9/;

    for my $interp_type ($gsl_interp_linear ,      #線形補間
                         $gsl_interp_cspline ){    #スプライン補間

        #この gsl_interp_accel_alloc() は十分に理解できていません
        my $accel = gsl_interp_accel_alloc();

        my $interp = gsl_interp_alloc($interp_type, scalar(@x));
        gsl_interp_init($interp, \@x, \@y, 10);

        for my $interp_x ( @interp_x_vals ){
            my $interp_y = gsl_interp_eval($interp,
                                           \@x, \@y,
                                           $interp_x,
                                           $accel);
            $interp_vals->{$interp_type}->{$interp_x} = $interp_y;
        }

        gsl_interp_free($interp);
        gsl_interp_accel_free($accel);
    }



    print "INTERPOLATED VALS\n";
    print "         LINEAR   CSPLINE\n";
    for my $interp_x ( @interp_x_vals ){
        print "X=", sprintf("%.2f",$interp_x);

        for my $interp_type ($gsl_interp_linear ,      #線形補間
                             $gsl_interp_cspline ){    #スプライン補間
            print " , Y=";
            print sprintf("%.2f",$interp_vals->{$interp_type}->{$interp_x});
        }
        print "\n";
    }
}

1;
__END__

↑と書くと、↓のように補間できます

$ ./foo.pl 
GENERATED ORIGINAL CURVE by sin() & cos()
X=0.00 , Y=1.00
X=1.42 , Y=1.54
X=2.45 , Y=1.35
X=3.07 , Y=2.09
X=3.62 , Y=3.04
X=4.52 , Y=5.99
X=5.86 , Y=5.87
X=7.33 , Y=7.30
X=8.49 , Y=8.39
X=9.21 , Y=9.78
INTERPOLATED VALS
         LINEAR   CSPLINE
X=0.00 , Y=1.00 , Y=1.00
X=1.00 , Y=1.38 , Y=1.54
X=1.50 , Y=1.53 , Y=1.52
X=3.50 , Y=2.83 , Y=2.76
X=5.00 , Y=5.95 , Y=6.40
X=7.20 , Y=7.17 , Y=7.15
X=8.60 , Y=8.60 , Y=8.56
X=9.00 , Y=9.37 , Y=9.33