Bspline derivative functions #40

Open
leto opened this Issue Oct 21, 2011 · 1 comment

Projects

None yet

1 participant

@leto
Owner
leto commented Oct 21, 2011

Feature request from Juergen Rose:

I am missing the gsl_bspline_deriv_workspace, gsl_bspline_derive_free,
gsl_bspline_deriv_eval and gsl_bspline_deriv_eval_nonzero functions in
Math::GSL. There are available at least in gsl-1.15. How are the chances
that they will be implemented in Math::GSL?

@leto
Owner
leto commented Oct 27, 2011
#!/usr/bin/perl -w

use strict;

use Math::GSL::BSpline qw/:all/;
use Math::GSL::Vector qw/:all/;
use Math::GSL::Matrix qw/:all/;
use Math::GSL::RNG qw/:all/;
use Math::GSL::Multifit qw /:all/;
use Math::GSL::Randist qw/:all/;
use Math::GSL::Statistics qw /:all/;

#  number of data points to fit
my $N=200;

#  number of fit coefficients
my $NCOEFFS=12;

#  nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4
my $NBREAK=$NCOEFFS -2;

my $n = $N;
my $ncoeffs = $NCOEFFS;
my $nbreak = $NBREAK;

my ($chisq, $Rsq, $dof, $tss);

my $seed=5;
#my $r = Math::GSL::RNG->new;
my $r = Math::GSL::RNG->new($gsl_rng_default,$seed);

#  allocate a cubic bspline workspace (k = 4)
my $bw = gsl_bspline_alloc(4, $nbreak);

my $B = Math::GSL::Vector->new($ncoeffs);
my $x = Math::GSL::Vector->new($n);
my $y = Math::GSL::Vector->new($n);
my $X = Math::GSL::Matrix->new($n, $ncoeffs);
my $c = Math::GSL::Vector->new($ncoeffs);
my $w = Math::GSL::Vector->new($n);
my $cov = Math::GSL::Matrix->new($ncoeffs, $ncoeffs);
my $mw = gsl_multifit_linear_alloc($n, $ncoeffs);

printf("#m=0,S=0\n");
#  this is the data to be fitted
for (my $i = 0; $i < $n; ++$i) {
    my $xi = (15.0 / ($N - 1)) * $i;
    my $yi = cos($xi) * exp(-0.1 * $xi);

    my $sigma = 0.1 * $yi;
    my $dy = gsl_ran_gaussian($r->raw(), $sigma);
    $yi += $dy;

    $x->set([$i], [$xi]);
    $y->set([$i], [$yi]);
    $w->set([$i], [1.0 / ($sigma * $sigma)]);

    printf("%f %f\n", $xi, $yi);
}

#  use uniform breakpoints on [0, 15]
gsl_bspline_knots_uniform(0.0, 15.0, $bw);

#  construct the fit matrix X
for (my $i = 0; $i < $n; ++$i) {
    my ($xi) = $x->get([$i]);

    #  compute B_j(xi) for all j
    gsl_bspline_eval($xi, $B->raw(), $bw);

    #  fill in row i of X
    for (my $j = 0; $j < $ncoeffs; ++$j) {
        my ($Bj) = $B->get([$j]);
        gsl_matrix_set($X->raw(), $i, $j, $Bj);
    }
}

#  do the fit
$chisq=gsl_multifit_wlinear($X->raw(), $w->raw(), $y->raw(), $c->raw(), $cov->raw(), $mw);

$dof = $n - $ncoeffs;

my @w_data=$w->as_list();
my @y_data=$y->as_list();
my $y_size=$y->length();
$tss = gsl_stats_wtss(\@w_data, 1, \@y_data, 1, $y->length());
$Rsq = 1.0 - $chisq / $tss;

printf STDERR "chisq/dof = %e, Rsq = %f\n",$chisq / $dof, $Rsq;

#  output the smoothed curve
{
    printf("#m=1,S=0\n");
    for (my $xi = 0.0; $xi < 15.0; $xi += 0.1) {
        gsl_bspline_eval($xi, $B->raw(), $bw);
        my ($rc,$yi,$yerr)=gsl_multifit_linear_est($B->raw(),$c->raw(),$cov->raw());
        printf("%f %f\n", $xi, $yi);
    }
}

my $dbw = gsl_bspline_deriv_alloc(4);

## Initialize dB with poison to ensure we overwrite it
my $dB = Math::GSL::Matrix->new($ncoeffs, 3);
#gsl_matrix_set_all($dB, GSL_NAN);

## output the smoothed curve via dB
{
    printf("#m=2,S=0\n");
    for (my $xi = 0.0; $xi < 15.0; $xi += 0.1) {
        gsl_bspline_deriv_eval($xi, 0, $dB, $bw, $dbw);
        ##  dB[i,j] = d^jB_i(x)/dx^j
        ##  vector dBk = get_column(dB,k)  k = k-th derivative ????
        my $dBk = gsl_matrix_column($dB,0);
        my ($rc,$yi,$yerr)=gsl_multifit_linear_est($dBk->raw(),$c->raw(),$cov->raw());
        printf("%f %f\n", $xi, $yi);
    }
}

## output the first derivative
{
    printf("#m=3,S=0\n");
    for (my $xi = 0.0; $xi < 15.0; $xi += 0.1) {
        gsl_bspline_deriv_eval($xi, 2, $dB, $bw, $dbw);
        ##  dB[i,j] = d^jB_i(x)/dx^j
        ##  vector dBk = get_column(dB,k)  k = k-th derivative ????
        my $dBk = gsl_matrix_column($dB,1);
        my ($rc,$yi,$yerr)=gsl_multifit_linear_est($dBk->raw(),$c->raw(),$cov->raw());
        printf("%f %f\n", $xi, $yi);
    }
}


$r->free();
gsl_bspline_deriv_free($dbw);
gsl_bspline_free($bw);
gsl_matrix_free($dB);
gsl_vector_free($B->raw());
gsl_vector_free($x->raw());
gsl_vector_free($y->raw());
gsl_matrix_free($X->raw());
gsl_vector_free($c->raw());
gsl_vector_free($w->raw());
gsl_matrix_free($cov->raw());
gsl_multifit_linear_free($mw);


#  The output can be plotted with gnu graph.
#
#    $ ./a.out > bspline.dat
#    chisq/dof = 1.118217e+00, Rsq = 0.989771
#    $ graph -T ps -X x -Y y -x 0 15 -y -1 1.3 < bspline.dat > bspline.ps
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment