gslでBスプラインでデータの平滑化

c

gslでBスプラインでデータの平滑化.
Bスプイラインは自然スプラインより計算が安定。
Bスプライン（拡張スプライン）同じところにノットを重ねることで
不連続を表現できる、３次だと、１回重ねると２回微分、２回重ねると１回微分、
３回重ねると関数そのものが不連続になる。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>

/* number of data points to fit */
#define N        1000

/* number of fit coefficients */
#define NCOEFFS 7

/* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
#define NBREAK   (NCOEFFS - 2)

void
spline(double* x_vector,double* y_vector, double* w_vector, int ndata,double* knots_vector,int nknots)
{
  const size_t n = (const size_t) ndata;
  const size_t ncoeffs = (const size_t) nknots+2;
  const size_t nbreak = (const size_t) nknots;
  size_t i, j;
  gsl_bspline_workspace *bw;
  gsl_vector *B;
  double dy;
  gsl_vector *c, *w;
  gsl_vector *x, *y;
  gsl_matrix *X, *cov;
  gsl_multifit_linear_workspace *mw;
  double chisq, Rsq, dof, tss;
 char* out_original_filename="out_original_data.dat";
 char* out_spline_filename="out_spline_data.dat";
  FILE* out_original;
  FILE* out_spline;

 // gsl_rng_env_setup();
//  r = gsl_rng_alloc(gsl_rng_default);

  /* allocate a cubic bspline workspace (k = 4) */
  bw = gsl_bspline_alloc(4, nbreak);
  B = gsl_vector_alloc(ncoeffs);

  x = gsl_vector_alloc(n);
  y = gsl_vector_alloc(n);
  X = gsl_matrix_alloc(n, ncoeffs);
  c = gsl_vector_alloc(ncoeffs);
  w = gsl_vector_alloc(n);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  mw = gsl_multifit_linear_alloc(n, ncoeffs);

 out_original=fopen(out_original_filename,"w");
 out_spline=fopen(out_spline_filename,"w");

// double *knots_vector;
// knots_vector=(double*)malloc(sizeof(double)*ncoeffs+2);

// double *x_vector;
// x_vector=(double*)malloc(sizeof(double)*n);

// double *y_vector;
// y_vector=(double*)malloc(sizeof(double)*n);

 //double *w_vector;
 // w_vector=(double*)malloc(sizeof(double)*n);


  gsl_vector *knots_vector2;
//input the data
/*
 for(i=0;i<n;i++){


      double sigma;
      double xi = i;
      double yi = 0.01*(i+1);

      sigma = 0.1 * yi;
      dy = gsl_ran_gaussian(r, sigma);
      yi += dy;
      y_vector[i]=yi;
      x_vector[i]=xi;
      w_vector[i]=1.0/sigma*sigma;
 }
*/
//Set the knots
/*
 for(i=0;i<ncoeffs-2;i++){
     double dx=(x_vector[n-1]+0.0)/((ncoeffs-2)-1.0);
     knots_vector[i]=i*dx;
     printf("knots %lf %lf\n",x_vector[n-1],knots_vector[i]);
 }
*/


  printf("#m=0,S=0\n");
  /* this is the data to be fitted */
  for (i = 0; i < n; ++i)
    {
      double wi=y_vector[i];
      double xi = x_vector[i];
      double yi = y_vector[i];

     // sigma = 0.1 * yi;
      //dy = gsl_ran_gaussian(r, sigma);
      //yi += dy;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(w, i, wi);

      fprintf(out_original,"%f %f\n", xi, yi);
    }
    knots_vector2=gsl_vector_alloc(ncoeffs-2);
    // gsl_vector_set(knots_vector2,0,min_x);
    for(i=0;i<ncoeffs-2;i++){
       gsl_vector_set(knots_vector2,i,knots_vector[i]);

     printf("%d\n",i);
     }
     // gsl_vector_set(knots_vector2,i,knots_vector[ncoeffs+1]);
  /* use uniform breakpoints on [0, 15] */
//  gsl_bspline_knots_uniform(0.0, gsl_vector_max(x), bw);
 printf("nbreak%d\n",bw->nbreak);
  gsl_bspline_knots(knots_vector2, bw);
  /* construct the fit matrix X */
  for (i = 0; i < n; ++i)
    {
      double xi = gsl_vector_get(x, i);

      /* compute B_j(xi) for all j */
      gsl_bspline_eval(xi, B, bw);

      /* fill in row i of X */
      for (j = 0; j < ncoeffs; ++j)
        {
          double Bj = gsl_vector_get(B, j);
          gsl_matrix_set(X, i, j, Bj);
        }
    }

  /* do the fit */
  gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);

  dof = n - ncoeffs;
  tss = gsl_stats_wtss(w->data, 1, y->data, 1, y->size);
  Rsq = 1.0 - chisq / tss;

  fprintf(stderr, "chisq/dof = %e, Rsq = %f\n",
                   chisq / dof, Rsq);

  /* output the smoothed curve */
  {
    double xi, yi, yerr;

    printf("#m=1,S=0\n");
    double delta=1.0;
    for (xi = 0.0; xi < gsl_vector_max(x); xi += delta)
      {
        gsl_bspline_eval(xi, B, bw);
        gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
        fprintf(out_spline,"%f %f\n", xi, yi);
      }
  }

//  gsl_rng_free(r);
  gsl_bspline_free(bw);
  gsl_vector_free(B);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_matrix_free(X);
  gsl_vector_free(c);
  gsl_vector_free(w);
  gsl_matrix_free(cov);
  gsl_multifit_linear_free(mw);
// free(knots_vector);
// free(x_vector);
// free(y_vector);
// free(w_vector);

 fclose(out_original);
 fclose(out_spline);

//  return 0;
}

int main(void){
        int i;
        int n=1000;
        int nknots=9;
        double dx,dy;

        //Initilaization of random variable
        gsl_rng *r;
        gsl_rng_env_setup();
        r = gsl_rng_alloc(gsl_rng_default);


        double *x_vector;
        x_vector=(double*)malloc(sizeof(double)*n);

        double *y_vector;
        y_vector=(double*)malloc(sizeof(double)*n);

        double *w_vector;
        w_vector=(double*)malloc(sizeof(double)*n);


        double *knots_vector;
        knots_vector=(double*)malloc(sizeof(double)*nknots);

        //Set the data
        for(i=0;i<n;i++){


             double sigma;
             double xi = i;
             double yi = 0.01*(i+1);

             sigma = 0.1 * yi;
             dy = gsl_ran_gaussian(r, sigma);
             yi += dy;
             y_vector[i]=yi;
             x_vector[i]=xi;
             w_vector[i]=1.0/sigma*sigma;
         }

        //Set the knots
         for(i=0;i<nknots;i++){
             double dx=(x_vector[n-1]+0.0)/((nknots)-1.0);
             knots_vector[i]=i*dx;
             printf("knots %lf %lf\n",x_vector[n-1],knots_vector[i]);
         }

        spline(x_vector,y_vector,w_vector,n,knots_vector,nknots);

         free(knots_vector);
         free(x_vector);
         free(y_vector);
         free(w_vector);


        return 0;


}