#include #include #include #include #include #include #include #include struct data { double *t; double *y; size_t n; }; double gaussian(const double a, const double b, const double c, const double t) { const double z = (t - b) / c; return (a * exp(-0.5 * z * z)); } //const double a = 1; //energy //const double b = 2.5; //sigema,distance //const double m = 12; /* m */ //const double n = 6; /* n */ //const double r0 = 0; /* cutoff */ //const double B = 1; // Extra attractive int func_f (const gsl_vector * x, void *params, gsl_vector * f) { struct data *d = (struct data *) params; double a = gsl_vector_get(x, 0); double b = gsl_vector_get(x, 1); double m = gsl_vector_get(x, 2); double n = gsl_vector_get(x, 3); double r0 = gsl_vector_get(x, 4); double B = gsl_vector_get(x, 5); size_t i; for (i = 0; i < d->n; ++i) { double ti = d->t[i]; double yi = d->y[i]; double y = a*(pow(b/(ti-r0),m)-B*pow(b/(ti-r0),n));//gaussian(a, b, c, ti); gsl_vector_set(f, i, yi - y); } return GSL_SUCCESS; } int func_df (const gsl_vector * x, void *params, gsl_matrix * J) { struct data *d = (struct data *) params; double a = gsl_vector_get(x, 0); double b = gsl_vector_get(x, 1); double m = gsl_vector_get(x, 2); double n = gsl_vector_get(x, 3); double r0 = gsl_vector_get(x, 4); double B = gsl_vector_get(x, 5); size_t i; for (i = 0; i < d->n; ++i) { double ti = d->t[i]; gsl_matrix_set(J, i, 0, pow((-b/(r0 - ti)),m) - B*pow((-b/(r0 - ti)),n)); gsl_matrix_set(J, i, 1, -a*((m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti) - (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)) ); gsl_matrix_set(J, i, 2, a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m)); gsl_matrix_set(J, i, 3, -B*a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n)); gsl_matrix_set(J, i, 4, a*((b*m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (B*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti))); gsl_matrix_set(J, i, 5, -a*pow((-b/(r0 - ti)),n)); } return GSL_SUCCESS; } int func_fvv (const gsl_vector * x, const gsl_vector * v, void *params, gsl_vector * fvv) { struct data *d = (struct data *) params; double a = gsl_vector_get(x, 0); double b = gsl_vector_get(x, 1); double m = gsl_vector_get(x, 2); double n = gsl_vector_get(x, 3); double r0 = gsl_vector_get(x, 4); double B = gsl_vector_get(x, 5); double va = gsl_vector_get(v, 0); double vb = gsl_vector_get(v, 1); double vm = gsl_vector_get(v, 2); double vn = gsl_vector_get(v, 3); double vr0 = gsl_vector_get(v, 4); double vB = gsl_vector_get(v, 5); size_t i; for (i = 0; i < d->n; ++i) { double ti = d->t[i]; //double Daa = 0; double Dab = (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti) - (m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti); double Dam = log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m); double Dan = -B*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n); double Dar0 = (b*m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (B*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti); double DaB = -pow((-b/(r0 - ti)),n); double Dbb = a*((m*pow((-b/(r0 - ti)),(m - 2))*(m - 1))/(r0 - ti)/(r0 - ti) - (B*n*pow((-b/(r0 - ti)),(n - 2))*(n - 1))/(r0 - ti)/(r0 - ti)); double Dbm = -a*(pow((-b/(r0 - ti)),(m - 1))/(r0 - ti) + (m*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)); double Dbn = a*((B*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti) + (B*n*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)); double Dbr0 = a*((m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti) - (b*m*pow((-b/(r0 - ti)),(m - 2))*(m - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti) + (B*b*n*pow((-b/(r0 - ti)),(n - 2))*(n - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti)); double DbB = (a*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti); double Dmm = a*log(-b/(r0 - ti))*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m); //double Dmn = 0; double Dmr0 = (a*b*m*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (a*pow((-b/(r0 - ti)),m))/(r0 - ti); //double DmB = 0; double Dnn = -B*a*log(-b/(r0 - ti))*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n); double Dnr0 = (B*a*pow((-b/(r0 - ti)),n))/(r0 - ti) - (B*a*b*n*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti); double DnB = -a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n); double Dr0r0 = -a*((2*b*m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti)/(r0 - ti) - (b*b*m*pow((-b/(r0 - ti)),(m - 2))*(m - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti)/(r0 - ti) - (2*B*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti)/(r0 - ti) + (B*b*b*n*pow((-b/(r0 - ti)),(n - 2))*(n - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti)/(r0 - ti)); double Dr0B = -(a*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti); //double DBB =0 double sum; sum = 2.0 * (va * vb * Dab + va * vm * Dam + va * vn * Dan + va * vr0 * Dar0 + va * vB * DaB) + 2.0 * (vb * vm * Dbm + vb * vn * Dbn + vb * vr0 * Dbr0 + vb * vB * DbB) + vb * vb * Dbb + vm * vm * Dmm + 2.0 * vm * vr0 * Dmr0 + vn * vn * Dnn + 2.0 * (vn * vr0 * Dnr0 + vn * vB * DnB + vr0 * vB * Dr0B) + vr0 * vr0 * Dr0r0; gsl_vector_set(fvv, i, sum); } return GSL_SUCCESS; } void callback(const size_t iter, void *params, const gsl_multifit_nlinear_workspace *w) { gsl_vector *f = gsl_multifit_nlinear_residual(w); gsl_vector *x = gsl_multifit_nlinear_position(w); double avratio = gsl_multifit_nlinear_avratio(w); double rcond; (void) params; /* not used */ /* compute reciprocal condition number of J(x) */ gsl_multifit_nlinear_rcond(&rcond, w); printf("calculating rcond, rcond=%f\n",rcond ); fprintf(stderr, "iter %2zu: a = %.4f, b = %.4f, m = %.4f ,n = %.4f, r0 = %.4f, B = %.4f, |a|/|v| = %.4f cond(J) = %8.4f, |f(x)| = %.4f\n", iter, gsl_vector_get(x, 0), gsl_vector_get(x, 1), gsl_vector_get(x, 2), gsl_vector_get(x, 3), gsl_vector_get(x, 4), gsl_vector_get(x, 5), avratio, 1.0 / rcond, gsl_blas_dnrm2(f)); } void solve_system(gsl_vector *x, gsl_multifit_nlinear_fdf *fdf, gsl_multifit_nlinear_parameters *params) { const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust; const size_t max_iter = 200; const double xtol = 1.0e-2; const double gtol = 1.0e-2; const double ftol = 1.0e-2; const size_t n = fdf->n; const size_t p = fdf->p; gsl_multifit_nlinear_workspace *work = gsl_multifit_nlinear_alloc(T, params, n, p); gsl_vector * f = gsl_multifit_nlinear_residual(work); gsl_vector * y = gsl_multifit_nlinear_position(work); int info; double chisq0, chisq, rcond; /* initialize solver */ gsl_multifit_nlinear_init(x, fdf, work); printf("ini done\n"); /* store initial cost */ gsl_blas_ddot(f, f, &chisq0); printf("store done\n"); /* iterate until convergence */ gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, callback, NULL, &info, work); printf("convergenced,info = %d\n" ,info); /* store final cost */ gsl_blas_ddot(f, f, &chisq); /* store cond(J(x)) */ gsl_multifit_nlinear_rcond(&rcond, work); gsl_vector_memcpy(x, y); /* print summary */ fprintf(stderr, "NITER = %zu\n", gsl_multifit_nlinear_niter(work)); fprintf(stderr, "NFEV = %zu\n", fdf->nevalf); fprintf(stderr, "NJEV = %zu\n", fdf->nevaldf); fprintf(stderr, "NAEV = %zu\n", fdf->nevalfvv); fprintf(stderr, "initial cost = %.12e\n", chisq0); fprintf(stderr, "final cost = %.12e\n", chisq); fprintf(stderr, "final x = (%.12e, %.12e, %12e)\n", gsl_vector_get(x, 0), gsl_vector_get(x, 1), gsl_vector_get(x, 2)); fprintf(stderr, "final cond(J) = %.12e\n", 1.0 / rcond); gsl_multifit_nlinear_free(work); } int main (void) { const size_t n = 300; /* number of data points to fit */ const size_t p = 6; /* number of model parameters */ const double a = 1; //energy const double b = 2.5; //sigema,distance const double m = 12; /* m */ const double n2 = 6; /* n */ const double r0 = 0; /* cutoff */ const double B = 1; // Extra attractive const gsl_rng_type * T = gsl_rng_default; gsl_vector *f = gsl_vector_alloc(n); gsl_vector *x = gsl_vector_alloc(p); gsl_multifit_nlinear_fdf fdf; gsl_multifit_nlinear_parameters fdf_params = gsl_multifit_nlinear_default_parameters(); struct data fit_data; gsl_rng * r; size_t i; gsl_rng_env_setup (); r = gsl_rng_alloc (T); fit_data.t = malloc(n * sizeof(double)); fit_data.y = malloc(n * sizeof(double)); fit_data.n = n; /* generate synthetic data with noise */ for (i = 0; i < n; ++i) { double t = (double)i / (double) n+1.5; double y0 = a*(pow(b/(t-r0),m)-B*pow(b/(t-r0),n2)); double dy = gsl_ran_gaussian (r, 0.1 * y0); fit_data.t[i] = t; fit_data.y[i] = y0 + dy*0; } /* define function to be minimized */ fdf.f = func_f; fdf.df = func_df; fdf.fvv = func_fvv; fdf.n = n; fdf.p = p; fdf.params = &fit_data; /* starting point */ gsl_vector_set(x, 0, 1.0); gsl_vector_set(x, 1, 2.4); gsl_vector_set(x, 2, 12.0); gsl_vector_set(x, 3, 6.0); gsl_vector_set(x, 4, 0.0); gsl_vector_set(x, 5, 1.0); fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel; solve_system(x, &fdf, &fdf_params); /* print data and model */ { double A = gsl_vector_get(x, 0); double B = gsl_vector_get(x, 1); double M = gsl_vector_get(x, 2); double N = gsl_vector_get(x, 3); double R0 = gsl_vector_get(x, 4); double BB = gsl_vector_get(x, 5); for (i = 0; i < n; ++i) { double ti = fit_data.t[i]; double yi = fit_data.y[i]; double fi = A*(pow(B/(ti-R0),M)-BB*pow(B/(ti-R0),N)); printf("%f %f %f\n", ti, yi, fi); } fprintf(stderr, "a = %.4f, b = %.4f, m = %.4f ,n = %.4f, r0 = %.4f, B = %.4f\n", A,//gsl_vector_get(x, 0), B,//gsl_vector_get(x, 1), M,//gsl_vector_get(x, 2), N,//gsl_vector_get(x, 3), R0,//gsl_vector_get(x, 4), BB/*gsl_vector_get(x, 5)*/); } gsl_vector_free(f); gsl_vector_free(x); gsl_rng_free(r); return 0; }