From 5032287c863cc770e7f720821716bdcd682dbcc9 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 17 Oct 2016 09:11:25 -0400
Subject: added some code to do data processing
---
makefile | 2 +-
src/anal_process.c | 119 ++++++++++++++++++++
src/fitting_functions.c | 290 ++++++++++++++++++++++++++++++++++++++++++++++++
3 files changed, 410 insertions(+), 1 deletion(-)
create mode 100644 src/anal_process.c
create mode 100644 src/fitting_functions.c
diff --git a/makefile b/makefile
index d0916ae..99440e9 100644
--- a/makefile
+++ b/makefile
@@ -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;
+}
+
--
cgit v1.2.3-70-g09d2