From 748552de6b57f6cfbe2f3668993f58657e241a57 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Fri, 8 Aug 2008 17:09:16 +0000 Subject: [PATCH] Fixed a few minor typos git-svn-id: svn://abax.physics.drexel.edu/sawsim/trunk@15 865a22a6-13cc-4084-8aa6-876098d8aa20 --- README | 6 +++- sawsim.nw | 95 +++++++++++++++++++++++++++++++++++++++---------------- 2 files changed, 73 insertions(+), 28 deletions(-) diff --git a/README b/README index 57aa974..32e9fb1 100644 --- a/README +++ b/README @@ -10,4 +10,8 @@ Run the unit tests with Depends on the GSL (GNU Scientific Library) development package. http://www.gnu.org/software/gsl/ -(libgsl0-dev in Debian-based distributions) + (libgsl0-dev in Debian-based distributions) + +The unit test programs check_* depend on the check package. + http://check.sourceforge.net/ + (check in Debian-based distributions) diff --git a/sawsim.nw b/sawsim.nw index 5a39047..133ac5c 100644 --- a/sawsim.nw +++ b/sawsim.nw @@ -827,7 +827,7 @@ int index_k_model(int n_models, k_model_getopt_t *models, char *name) &num_param_args, ¶m_args); \ if (num_param_args != model->num_params) { \ fprintf(stderr, \ - "%s model %s expected %d params," \ + "%s model %s expected %d params,\n", \ role, model->name, model->num_params); \ fprintf(stderr, \ "not the %d params in '%s'\n", \ @@ -1868,8 +1868,8 @@ double tension_balance(int num_tension_groups, { static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL}; double F, xo_guess, xo, lb, ub; - double min_dx=1e-15, min_dy=1e-16; - int max_steps=100, i; + double min_relx=1e-6, min_rely=1e-6; + int max_steps=200, i; assert(num_tension_groups > 0); assert(tension_handlers != NULL); @@ -1921,7 +1921,8 @@ double tension_balance(int num_tension_groups, fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n", __FUNCTION__, x, lb, ub); #endif - xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_dx, min_dy, max_steps, NULL); + xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL); + F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */ #ifdef DEBUG if (num_tension_groups > 1) { fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x); @@ -1954,42 +1955,60 @@ double x_of_xo(double xo, void *pdata) x_of_xo_data_t *data = (x_of_xo_data_t *)pdata; double F, x=0, xi, xi_guess, lb, ub; int i; - double min_dx=1e-15, min_dy=1e-16; - int max_steps=100; + double min_relx=1e-6, min_rely=1e-6; + int max_steps=200; assert(data->n_groups > 0); data->xi[0] = xo; F = (*data->tension_handler[0])(xo, data->handler_data[0]); +#ifdef DEBUG + fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F); +#endif x += xo; if (data->xi) data->xi[0] = xo; for (i=1; i < data->n_groups; i++) { if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */ else xi_guess = data->xi[i]; +#ifdef DEBUG + fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F); +#endif oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub); - xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_dx, min_dy, max_steps, NULL); + xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL); +#ifdef DEBUG + fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F); +#endif data->xi[i] = xi; x += xi; if (data->xi) data->xi[i] = xi; } +#ifdef DEBUG + fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x); +#endif return x; } @ -Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited number of steps for monotonically increasing $f(x)$. -Simple bisection, so it's robust and converges fairly quickly. +Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited +number of steps for monotonically increasing $f(x)$. Simple +bisection, so it's robust and converges fairly quickly. You can set +the maximum number of steps to take, but if you have enough processor +power, it's better to limit your precision with the relative limits +[[min_relx]] and [[y]]. These ensure that the width of the bracket is +small on the length scale of it's larger side. Note that \emph{both} +relative limits must be satisfied for the search to stop. <>= /* equivalent to gsl_func_t */ typedef double one_dim_fn_t(double x, void *params); @ <>= double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, - double min_dx, double min_dy, int max_steps, + double min_relx, double min_rely, int max_steps, double *pYx); @ <>= /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */ double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, - double min_dx, double min_dy, int max_steps, + double min_relx, double min_rely, int max_steps, double *pYx) { double yx, ylb, yub, x; @@ -2016,7 +2035,7 @@ double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, #ifdef DEBUG fprintf(stderr, "%s:\tbracket: lb %g, x %g, ub %g\tylb %g, yx %g, yub %g\t(y %g)\n", __FUNCTION__, lb, x, ub, ylb, yx, yub, y); #endif - if (ub-lb < min_dx || yub-ylb < min_dy || yx == y) break; + if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break; else if (yx > y) { ub=x; yub=yx; } else /* yx < y */ { lb=x; ylb=yx; } } @@ -2128,7 +2147,7 @@ Suite *test_suite (void) <>= <> -double hooke(void *pK, double x) +double hooke(double x, void *pK) { assert(x >= 0); return *((double*)pK) * x; @@ -2138,11 +2157,8 @@ START_TEST(test_single_function) { double k=5, x=3, last_x=2.0, F; one_dim_fn_t *handlers[] = {&hooke}; - void *data[] = {&k}; - double xi[] = {0.0}; - int active[] = {1}; - int new_active_groups = 1; - F = tension_balance(1, handlers, data, active, new_active_groups, xi, last_x, x); + void *data[] = {&k}; + F = tension_balance(1, handlers, data, last_x, x); fail_unless(F = k*x, NULL); } END_TEST @@ -2154,18 +2170,15 @@ The analytic solution for a general number of springs is given in Appendix \ref{ START_TEST(test_double_hooke) { double k1=5, k2=4, x=3, last_x=2.0, F, Fe, x1e, x2e; - one_dim_fn_t *handlers[] = {&hooke, &hooke, NULL}; - void *data[] = {&k1, &k2, NULL}; - double xi[] = {0.0, 0.0, 0.0}; - int active[] = {1, 1, 0}; - int new_active_groups = 1; - F = tension_balance(3, handlers, data, active, new_active_groups, xi, last_x, x); + one_dim_fn_t *handlers[] = {&hooke, &hooke}; + void *data[] = {&k1, &k2}; + F = tension_balance(2, handlers, data, last_x, x); x1e = x*k2/(k1+k2); Fe = k1*x1e; x2e = x1e*k1/k2; //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F); - CHECK_ERR(1e-6, x1e, xi[0]); - CHECK_ERR(1e-6, x2e, xi[1]); + //CHECK_ERR(1e-6, x1e, xi[0]); + //CHECK_ERR(1e-6, x2e, xi[1]); CHECK_ERR(1e-6, Fe, F); } END_TEST @@ -3149,6 +3162,9 @@ double const_tension_handler(double x, void *pdata) list = head(list); assert(list != NULL); /* empty active group?! */ F = ((const_tension_param_t *)list->d)->F; +#ifdef DEBUG + fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F); +#endif while (list != NULL) { assert(((const_tension_param_t *)list->d)->F == F); list = list->next; @@ -3248,9 +3264,17 @@ double hooke_handler(double x, void *pdata) while (list != NULL) { assert( ((hooke_param_t *)list->d)->k > 0 ); k += 1.0/ ((hooke_param_t *)list->d)->k; +#ifdef DEBUG + fprintf(stderr, "%s: Hooke spring %g N/m in series\n", + __FUNCTION__, ((hooke_param_t *)list->d)->k); +#endif list = list->next; } k = 1.0 / k; +#ifdef DEBUG + fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n", + __FUNCTION__, k, x, k*x); +#endif return k*x; } @ @@ -3364,8 +3388,16 @@ double wlc_handler(double x, void *pdata) /* enforce consistent p */ assert( ((wlc_param_t *)list->d)->p == p); L += ((wlc_param_t *)list->d)->L; +#ifdef DEBUG + fprintf(stderr, "%s: WLC section %g m long in series\n", + __FUNCTION__, ((wlc_param_t *)list->d)->L); +#endif list = list->next; } +#ifdef DEBUG + fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n", + __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L)); +#endif return wlc(x, data->env->T, p, L); } @ @@ -3497,8 +3529,16 @@ double fjc_handler(double x, void *pdata) /* enforce consistent l */ assert( ((fjc_param_t *)list->d)->l == l); N += ((fjc_param_t *)list->d)->N; +#ifdef DEBUG + fprintf(stderr, "%s: FJC section %d links long in series\n", + __FUNCTION__, ((fjc_param_t *)list->d)->N); +#endif list = list->next; } +#ifdef DEBUG + fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n", + __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N)); +#endif return fjc(x, data->env->T, l, N); } @ @@ -3556,6 +3596,7 @@ char fjc_param_string[]="0.5e-9,200"; \subsection{Tension model implementation} <>= +//#define DEBUG <> <> #include "tension_model.h" @@ -4243,7 +4284,7 @@ The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate It works so long as [[val_1]] is not [[false]]. <>= -{"kramers", "thermalized, diffusion-limited transitions (function parameters are parsed by Python. For distances, the environment variables force 'F' and temperature 'T' are defined, as well as the Boltzmann constant 'kB'. For energies, the position 'x' is also defined. Functions must evaluate to a float representing their output, and produce no output on stdout.", &kramers_k, num_kramers_params, kramers_param_descriptions, (char *)kramers_param_string, &string_create_kramers_param_t, &destroy_kramers_param_t} +{"kramers", "thermalized, diffusion-limited transitions (function parameters are parsed by Python. For distances, the environment variables force 'F' and temperature 'T' are defined, as well as the Boltzmann constant 'kB'. For energies, the position 'x' is also defined. Functions must evaluate to a float representing their output and produce no output on stdout.", &kramers_k, num_kramers_params, kramers_param_descriptions, (char *)kramers_param_string, &string_create_kramers_param_t, &destroy_kramers_param_t} @ \subsection{Kramer's model (natural cubic splines)} -- 2.26.2