diff options
| -rw-r--r-- | makefile | 2 | ||||
| -rw-r--r-- | src/anal_process.c | 119 | ||||
| -rw-r--r-- | src/fitting_functions.c | 290 |
3 files changed, 410 insertions, 1 deletions
@@ -4,7 +4,7 @@ CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-fi LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler #-ltcmalloc OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand -BIN = corr_test fracture +BIN = corr_test fracture fitting_functions anal_process all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} diff --git a/src/anal_process.c b/src/anal_process.c new file mode 100644 index 0000000..ed2f1f7 --- /dev/null +++ b/src/anal_process.c @@ -0,0 +1,119 @@ + +#include "fracture.h" +#include <gsl/gsl_sf_erf.h> +#include <gsl/gsl_sf_laguerre.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_multifit_nlinear.h> + +void mom(uint_t len, uint_t n, uint32_t *data, double result[2]) { + uint_t total = 0; + double mom = 0; + double mom_err = 0; + + for (uint_t i = 0; i < len; i++) { + uint32_t datai = data[i]; + double in = pow(i, n); + + total += datai; + mom += in * datai; + mom_err += gsl_pow_2(in) * datai; + } + + double momf = mom / total; + double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total); + + result[0] = momf; + result[1] = momf_err; +} + +int main(int argc, char *argv[]) { + uint_t nc = argc - 1; + uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); + double *betas_c = (double *)malloc(nc * sizeof(double)); + double *vals_c2 = (double *)malloc(nc * sizeof(double)); + double *errors_c2 = (double *)malloc(nc * sizeof(double)); + double *vals_c3 = (double *)malloc(nc * sizeof(double)); + double *errors_c3 = (double *)malloc(nc * sizeof(double)); + + char *out_filename = (char *)malloc(100 * sizeof(char)); + uint_t out_filename_len = 0; + + for (uint_t i = 0; i < nc; i++) { + char *fn = argv[1 + i]; + char *buff = (char *)malloc(20 * sizeof(char)); + uint_t pos = 0; uint_t c = 0; + while (c < 4) { + if (fn[pos] == '_') { + c++; + } + if (i == 0) { + out_filename[pos] = fn[pos]; + out_filename_len++; + } + pos++; + } + c = 0; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + Ls_c[i] = atoi(buff); + c = 0; + pos++; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + betas_c[i] = atof(buff); + free(buff); + + uint_t dist_len = 4 * pow(Ls_c[i], 2); + uint32_t *dist = malloc(dist_len * sizeof(uint32_t)); + FILE *dist_file = fopen(fn, "rb"); + fread(dist, sizeof(uint32_t), dist_len, dist_file); + fclose(dist_file); + { + double mom2[2]; + mom(dist_len, 2, dist, mom2); + vals_c2[i] = mom2[0]; + errors_c2[i] = mom2[1]; + + double mom3[2]; + mom(dist_len, 3, dist, mom3); + vals_c3[i] = mom3[0]; + errors_c3[i] = mom3[1]; + } + free(dist); + } + + out_filename[out_filename_len-1] = '.'; + out_filename[out_filename_len] = 't'; + out_filename[out_filename_len+1] = 'x'; + out_filename[out_filename_len+2] = 't'; + out_filename[out_filename_len+3] = '\0'; + + FILE *out_file = fopen(out_filename, "w"); + + for (uint_t i = 0; i < nc; i++) { + fprintf(out_file, "%u %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i]); + } + + fclose(out_file); + + + free(Ls_c); + free(betas_c); + free(vals_c2); + free(errors_c2); + free(vals_c3); + free(errors_c3); + + return 0; +} + diff --git a/src/fitting_functions.c b/src/fitting_functions.c new file mode 100644 index 0000000..aa7a305 --- /dev/null +++ b/src/fitting_functions.c @@ -0,0 +1,290 @@ + +#include "fracture.h" +#include <gsl/gsl_sf_erf.h> +#include <gsl/gsl_sf_laguerre.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_multifit_nlinear.h> + +struct model_params { + double nu; + double ac; + double mc; + double ac02; + double ac03; + double ac12; + double ac13; + double bc02; + double bc03; + double bc12; + double bc13; + double Ac20; + double Ac21; + double Ac22; + double Ac23; + double Ac30; + double Ac31; + double Ac32; + double Ac33; + double Bc20; + double Bc21; + double Bc22; + double Bc23; + double Bc30; + double Bc31; + double Bc32; + double Bc33; +}; + +struct data { + uint_t nc; + uint_t *Ls_c; + double *betas_c; + double *vals_c2; + double *vals_c3; +}; + +double Jc(uint_t n, uint_t o, struct model_params *par, double x) { + double nu, ac, mc, ac0n, bc0n, *Acn, yc, sum; + + nu = par->nu; + ac = par->ac; + mc = par->mc; + ac0n = *(&(par->ac02) + (n - 2) * sizeof(double)); + bc0n = *(&(par->bc02) + (n - 2) * sizeof(double)); + Acn = &(par->Ac20) + (n - 2) * o * sizeof(double); + + yc = (log(x) - mc) / ac; + + sum = 0; + + for (uint_t i = 0; i < o; i++) { + sum += Acn[i] * gsl_sf_laguerre_n(i, 0, yc); + } + + return exp(ac0n + bc0n * (2 - gsl_sf_erf(yc)) + exp(-gsl_pow_2(yc)) * sum); +} + +double Kc(uint_t n, uint_t o, struct model_params *par, double x) { + double nu, ac, mc, ac1n, bc1n, *Bcn, yc, sum; + + nu = par->nu; + ac = par->ac; + mc = par->mc; + ac1n = *(&(par->ac12) + (n - 2) * sizeof(double)); + bc1n = *(&(par->bc12) + (n - 2) * sizeof(double)); + Bcn = &(par->Bc20) + (n - 2) * o * sizeof(double); + + yc = (log(x) - mc) / ac; + + sum = 0; + + for (uint_t i = 0; i < o; i++) { + sum += Bcn[i] * gsl_sf_laguerre_n(i, 0, yc); + } + + return exp(ac1n + bc1n * (2 - gsl_sf_erf(yc)) + exp(-gsl_pow_2(yc)) * sum); +} + +double sc(uint_t n, uint_t o, struct model_params *par, uint_t L, double beta) { + double nu, bLnu, deltanu, sigmanu, tau, Lntau, Ldelta; + + nu = par->nu; + deltanu = 1.5; + sigmanu = 48. / 91; + tau = 187. / 91; + + bLnu = beta * exp(log(L) / nu); + Lntau = exp((n + 1 - tau) / sigmanu * log(L)); + Ldelta = exp(-deltanu * log(L)); + + return Lntau * (Jc(n, o, par, bLnu) + Ldelta * Kc(n, 0, par, bLnu)); +} + +int f_moms(const gsl_vector *x, void *params, gsl_vector *f) { + struct data *dat = (struct data *)params; + + for (uint_t i = 0; i < dat->nc; i++) { + double f2i = sc(2, 4, (struct model_params *)x->data, dat->Ls_c[i], dat->betas_c[i]); + double F2i = dat->vals_c2[i]; + + f2i = log(f2i); + F2i = log(F2i); + + gsl_vector_set(f, i, f2i - F2i); + + double f3i = sc(3, 4, (struct model_params *)x->data, dat->Ls_c[i], dat->betas_c[i]); + double F3i = dat->vals_c3[i]; + + f3i = log(f3i); + F3i = log(F3i); + + gsl_vector_set(f, dat->nc + i, f3i - F3i); + } + + return GSL_SUCCESS; +} + +void mom(uint_t len, uint_t n, uint32_t *data, double result[2]) { + uint_t total = 0; + double mom = 0; + double mom_err = 0; + + for (uint_t i = 0; i < len; i++) { + uint32_t datai = data[i]; + double in = pow(i, n); + + total += datai; + mom += in * datai; + mom_err += gsl_pow_2(in) * datai; + } + + double momf = mom / total; + double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total); + + result[0] = momf; + result[1] = momf_err; +} + +void +callback(const size_t iter, void *params, + const gsl_multifit_nlinear_workspace *w) +{ + gsl_vector *f = gsl_multifit_nlinear_residual(w); + + fprintf(stderr, "iter %2zu: |f(x)| = %.4f\n", + iter, + gsl_blas_dnrm2(f)); +} + + +int main(int argc, char *argv[]) { + struct data *d = (struct data *)malloc(sizeof(struct data)); + uint_t nc = argc - 1; + uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); + double *betas_c = (double *)malloc(nc * sizeof(double)); + double *vals_c2 = (double *)malloc(nc * sizeof(double)); + double *weights_c2 = (double *)malloc(nc * sizeof(double)); + double *vals_c3 = (double *)malloc(nc * sizeof(double)); + double *weights_c3 = (double *)malloc(nc * sizeof(double)); + + double chisq, chisq0; + + d->nc = nc; + d->Ls_c = Ls_c; + d->betas_c = betas_c; + d->vals_c2 = vals_c2; + d->vals_c3 = vals_c3; + + for (uint_t i = 0; i < nc; i++) { + char *fn = argv[1 + i]; + char *buff = malloc(20 * sizeof(char)); + uint_t pos = 0; uint_t c = 0; + while (c < 4) { + if (fn[pos] == '_') { + c++; + } + pos++; + } + c = 0; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + Ls_c[i] = atoi(buff); + c = 0; + pos++; + while (fn[pos] != '_') { + buff[c] = fn[pos]; + pos++; + c++; + } + buff[c] = '\0'; + betas_c[i] = atof(buff); + free(buff); + + uint_t dist_len = 4 * pow(Ls_c[i], 2); + uint32_t *dist = malloc(dist_len * sizeof(uint32_t)); + FILE *dist_file = fopen(fn, "rb"); + fread(dist, sizeof(uint32_t), dist_len, dist_file); + fclose(dist_file); + { + double mom2[2]; + mom(dist_len, 2, dist, mom2); + vals_c2[i] = mom2[0]; + weights_c2[i] = mom2[1]; + + double mom3[2]; + mom(dist_len, 3, dist, mom3); + vals_c3[i] = mom3[0]; + weights_c3[i] = mom3[1]; + } + free(dist); + } + + const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust; + gsl_multifit_nlinear_workspace *w; + gsl_multifit_nlinear_fdf fdf; + gsl_multifit_nlinear_parameters fdf_params = gsl_multifit_nlinear_default_parameters(); + uint_t n = 2 * nc; + uint_t p = 27; + + double x_init[27] = { 1.56, .3, 2, -6, -10, -10, -10, 10, 15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + double weights[n]; + for (uint_t i = 0; i < nc; i++) { + weights[i] = 1/gsl_pow_2(weights_c2[i] / vals_c2[i]); + weights[nc + i] = 1/gsl_pow_2(weights_c3[i] / vals_c3[i]); + } + + gsl_vector_view x = gsl_vector_view_array(x_init, p); + gsl_vector_view wts = gsl_vector_view_array(weights, n); + + gsl_vector *f; + + fdf.f = f_moms; + fdf.df = NULL; + fdf.fvv = NULL; + fdf.n = n; + fdf.p = p; + fdf.params = d; + + fdf_params.trs = gsl_multifit_nlinear_trs_lm; + + w = gsl_multifit_nlinear_alloc(T, &fdf_params, n, p); + gsl_multifit_nlinear_winit(&x.vector, &wts.vector, &fdf, w); + f = gsl_multifit_nlinear_residual(w); + gsl_blas_ddot(f, f, &chisq0); + + const double xtol = 0.0; + const double gtol = 0.0; + const double ftol = 0.0; + + int info; + int status = gsl_multifit_nlinear_driver(100, xtol, gtol, ftol, callback, NULL, &info, w); + gsl_blas_ddot(f, f, &chisq); + + fprintf(stderr, "summary from method '%s/%s'\n", + gsl_multifit_nlinear_name(w), + gsl_multifit_nlinear_trs_name(w)); + + fprintf(stderr, "number of iterations: %zu\n", + gsl_multifit_nlinear_niter(w)); + fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf); + fprintf(stderr, "reason for stopping: %s\n", + (info == 1) ? "small step size" : "small gradient"); + fprintf(stderr, "initial |f(x)| = %g\n", sqrt(chisq0)); + fprintf(stderr, "final |f(x)| = %g\n", sqrt(chisq)); + + free(Ls_c); + free(betas_c); + free(vals_c2); + free(weights_c2); + free(vals_c3); + free(weights_c3); + + return 0; +} + |