histogram.c in gsl-2.1.0

- old
+ new

@@ -1167,10 +1167,13 @@
   int status;
   size_t iter = 0, binstart, binend;
   size_t n, dof;      /* # of data points */
   size_t p = 3;  /* # of fitting parameters */
   gsl_multifit_function_fdf f;
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  gsl_matrix *J = NULL;
+#endif
   gsl_matrix *covar = NULL;
   gsl_vector *x = NULL;
   double sigma, mean, height, errs, errm, errh, chi2;
   Data_Get_Struct(obj, gsl_histogram, h);
   binstart = 0;
@@ -1195,10 +1198,13 @@
   hh.h = h;
   hh.binstart = binstart;
   hh.binend = binend;
   n = binend - binstart + 1;
 
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  J = gsl_matrix_alloc(n, p);
+#endif
   covar = gsl_matrix_alloc(p, p);
 
   f.f = Gaussian_f;
   f.df = Gaussian_df;
   f.fdf = Gaussian_fdf;
@@ -1217,19 +1223,27 @@
     status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
   } while (status == GSL_CONTINUE);
   sigma = sqrt(gsl_vector_get(s->x, 0));
   mean = gsl_vector_get(s->x, 1);
   height = gsl_vector_get(s->x, 2)*sigma*sqrt(2*M_PI);
+#ifdef GSL_MULTIFIT_FDFSOLVER_J
   gsl_multifit_covar(s->J, 0.0, covar);
+#else
+  gsl_multifit_fdfsolver_jac(s, J);
+  gsl_multifit_covar(J, 0.0, covar);
+#endif
   chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));   /* not reduced chi-square */
   dof = n - p;
   errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0))/sigma/2;
   errm = sqrt(chi2/dof*gsl_matrix_get(covar, 1, 1));
   errh = sqrt(chi2/dof*gsl_matrix_get(covar, 2, 2));
 
   gsl_multifit_fdfsolver_free(s);
   gsl_vector_free(x);
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  gsl_matrix_free(J);
+#endif
   gsl_matrix_free(covar);
   return rb_ary_new3(8, rb_float_new(sigma), rb_float_new(mean),
                      rb_float_new(height), rb_float_new(errs),
                      rb_float_new(errm), rb_float_new(errh),
                      rb_float_new(chi2), INT2FIX(dof));
@@ -1303,10 +1317,13 @@
   int status;
   size_t iter = 0, binstart, binend;
   size_t n, dof;      /* # of data points */
   size_t p = 2;  /* # of fitting parameters */
   gsl_multifit_function_fdf f;
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  gsl_matrix *J = NULL;
+#endif
   gsl_matrix *covar = NULL;
   gsl_vector *x = NULL;
   double sigma, height, errs, errh, chi2;
   Data_Get_Struct(obj, gsl_histogram, h);
   binstart = 0;
@@ -1330,10 +1347,13 @@
   hh.h = h;
   hh.binstart = binstart;
   hh.binend = binend;
   n = binend - binstart + 1;
 
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  J = gsl_matrix_alloc(n, p);
+#endif
   covar = gsl_matrix_alloc(p, p);
 
   f.f = Rayleigh_f;
   f.df = Rayleigh_df;
   f.fdf = Rayleigh_fdf;
@@ -1351,18 +1371,26 @@
     if (status) break;
     status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
   } while (status == GSL_CONTINUE);
   sigma = sqrt(gsl_vector_get(s->x, 0));
   height = gsl_vector_get(s->x, 1)*sigma*sigma;
+#ifdef GSL_MULTIFIT_FDFSOLVER_J
   gsl_multifit_covar(s->J, 0.0, covar);
+#else
+  gsl_multifit_fdfsolver_jac(s, J);
+  gsl_multifit_covar(J, 0.0, covar);
+#endif
   chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));   /* not reduced chi-square */
   dof = n - p;
   errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0))/sigma/2;
   errh = sqrt(chi2/dof*gsl_matrix_get(covar, 1, 1));
 
   gsl_multifit_fdfsolver_free(s);
   gsl_vector_free(x);
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  gsl_matrix_free(J);
+#endif
   gsl_matrix_free(covar);
   return rb_ary_new3(6, rb_float_new(sigma),
                      rb_float_new(height), rb_float_new(errs),
                      rb_float_new(errh),
                      rb_float_new(chi2), INT2FIX(dof));
@@ -1439,10 +1467,13 @@
   int status;
   size_t iter = 0, binstart, binend;
   size_t n, dof;      /* # of data points */
   size_t p = 2;  /* # of fitting parameters */
   gsl_multifit_function_fdf f;
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  gsl_matrix *J = NULL;
+#endif
   gsl_matrix *covar = NULL;
   gsl_vector *x = NULL;
   double b, height, errs, errh, chi2;
   Data_Get_Struct(obj, gsl_histogram, h);
   binstart = 0;
@@ -1466,10 +1497,13 @@
   hh.h = h;
   hh.binstart = binstart;
   hh.binend = binend;
   n = binend - binstart + 1;
 
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  J = gsl_matrix_alloc(n, p);
+#endif
   covar = gsl_matrix_alloc(p, p);
 
   f.f = xExponential_f;
   f.df = xExponential_df;
   f.fdf = xExponential_fdf;
@@ -1487,17 +1521,25 @@
     if (status) break;
     status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
   } while (status == GSL_CONTINUE);
   b = gsl_vector_get(s->x, 0);
   height = gsl_vector_get(s->x, 1);
+#ifdef GSL_MULTIFIT_FDFSOLVER_J
   gsl_multifit_covar(s->J, 0.0, covar);
+#else
+  gsl_multifit_fdfsolver_jac(s, J);
+  gsl_multifit_covar(J, 0.0, covar);
+#endif
   chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));   /* not reduced chi-square */
   dof = n - p;
   errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0));
   errh = sqrt(chi2/dof*gsl_matrix_get(covar, 1, 1));
 
   gsl_multifit_fdfsolver_free(s);
   gsl_vector_free(x);
+#ifndef GSL_MULTIFIT_FDFSOLVER_J
+  gsl_matrix_free(J);
+#endif
   gsl_matrix_free(covar);
   return rb_ary_new3(6, rb_float_new(b),
                      rb_float_new(height), rb_float_new(errs),
                      rb_float_new(errh),
                      rb_float_new(chi2), INT2FIX(dof));