Class: MultiFit::FdfSolver

# File 'ext/multifit.c'

/* Singleton method */
static VALUE rb_gsl_multifit_fit(int argc, VALUE *argv, VALUE klass)
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *solver;
  int status;
  size_t iter = 0, i;
  size_t n, dof;      /* # of data points */
  size_t p;           /* # of fitting parameters */
  gsl_multifit_function_fdf f;
  gsl_matrix *covar = NULL;
  gsl_vector *v = NULL;
  gsl_vector *x, *y, *w = NULL;
  gsl_vector_view xx, yy, ww;
  gsl_vector *vout, *verr;
  int flag = 0;
  double chi2;
  char fittype[256];
  struct fitting_xydata xydata;
  if (argc < 3) rb_raise(rb_eArgError, "too few arguments");
  switch (TYPE(argv[argc-1])) {
  case T_ARRAY:
    v = get_vector(argv[argc-1]);
    flag = 1;
    argc--;
    break;
  case T_STRING:
    /* do nothing */
    break;
  default:
    Data_Get_Vector(argv[argc-1], v);
    flag = 0;
    argc--;
    break;
  }
  x = &xx.vector;
  y = &yy.vector;
  w = &ww.vector;
  switch (argc) {
  case 3:
    Data_Get_Vector(argv[0], x);
    Data_Get_Vector(argv[1], y);
    w = NULL;
    strcpy(fittype, STR2CSTR(argv[2]));
    break;
  case 4:
    Data_Get_Vector(argv[0], x);
    Data_Get_Vector(argv[1], w);
    Data_Get_Vector(argv[2], y);
    strcpy(fittype, STR2CSTR(argv[3]));
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments");
    break;
  }

  xydata.x = x;
  xydata.y = y;
  xydata.w = w;
  n = x->size;

  set_fittype(&f, fittype, &p, &v, &flag);

  f.n = n;
  f.p = p;
  f.params = &xydata;

  T = gsl_multifit_fdfsolver_lmsder;
  solver = gsl_multifit_fdfsolver_alloc(T, n, p);
  gsl_multifit_fdfsolver_set(solver, &f, v);

  do {
    iter++;
    status = gsl_multifit_fdfsolver_iterate(solver);
    if (status) break;
    status = gsl_multifit_test_delta(solver->dx, solver->x, 1e-6, 1e-6);
  } while (status == GSL_CONTINUE);
  vout = gsl_vector_alloc(p);
  verr = gsl_vector_alloc(p);
  gsl_vector_memcpy(vout, solver->x);
  covar = gsl_matrix_alloc(p, p);
  chi2 = gsl_pow_2(gsl_blas_dnrm2(solver->f));   /* not reduced chi-square */
  dof = n - p;
  gsl_multifit_covar(solver->J, 0.0, covar);
  for (i = 0; i < p; i++)
    gsl_vector_set(verr, i, sqrt(chi2/dof*gsl_matrix_get(covar, i, i)));
  gsl_matrix_free(covar);
  if (flag == 1) gsl_vector_free(v);
  gsl_multifit_fdfsolver_free(solver);
  return rb_ary_new3(4,
             Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, vout),
             Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, verr),
             rb_float_new(chi2), INT2FIX(dof));
}

# File 'ext/multifit.c'

/**********/

static VALUE rb_gsl_multifit_fdfsolver_new(int argc, VALUE *argv, VALUE klass)
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *solver = NULL;
  size_t n, p;
  char name[64];
  switch (argc) {
  case 3:
    switch (TYPE(argv[0])) {
    case T_STRING:
      strcpy(name, STR2CSTR(argv[0]));
      if (str_tail_grep(name, "lmsder") == 0) {
    T = gsl_multifit_fdfsolver_lmsder;
      } else if (str_tail_grep(name, "lmder") == 0) {
    T = gsl_multifit_fdfsolver_lmder;
      } else {
    rb_raise(rb_eTypeError, "unknown solver type %s (lmsder of lmder)",
         name);
      }
      break;
    case T_FIXNUM:
      switch (FIX2INT(argv[0])) {
      case GSL_MULTIFIT_FDFSOLVER_LMSDER:
    T = gsl_multifit_fdfsolver_lmsder;
    break;
      case GSL_MULTIFIT_FDFSOLVER_LMDER:
    T = gsl_multifit_fdfsolver_lmder;
    break;
      default:
    rb_raise(rb_eTypeError, 
         "unknown solver type (GSL::MultiFit::FdfSolver::LMSDER or LMDER expected)");
    break;
      }
      break;
    default:
      rb_raise(rb_eTypeError, "wrong argument type %s (Fixnum or String)",
           rb_class2name(CLASS_OF(argv[0])));
      break;
    }
    CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
    n = FIX2INT(argv[1]);
    p = FIX2INT(argv[2]);
    break;
  case 2:
    CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]);
    T = gsl_multifit_fdfsolver_lmsder;
    n = FIX2INT(argv[0]);
    p = FIX2INT(argv[1]);
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments");
    break;
  }
  solver = gsl_multifit_fdfsolver_alloc(T, n, p);
  return Data_Wrap_Struct(klass, 0, gsl_multifit_fdfsolver_free, solver);
}