From 0dfa0698ae539317832964d6e785afd7b3242361 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Tue, 24 Mar 2009 13:40:02 -0400 Subject: [PATCH] Added inverse tension functions for speed. Bumped to version 0.7. --- src/sawsim.bib | 10 ++ src/sawsim.nw | 298 +++++++++++++++++++++++++------ testing/bell_rate/run_test.sh | 27 ++- testing/common/fit.pyc | Bin 4361 -> 0 bytes testing/const_rate/const_rate.sh | 20 ++- 5 files changed, 284 insertions(+), 71 deletions(-) delete mode 100644 testing/common/fit.pyc diff --git a/src/sawsim.bib b/src/sawsim.bib index 142ab7c..6977cdf 100644 --- a/src/sawsim.bib +++ b/src/sawsim.bib @@ -3807,3 +3807,13 @@ doi = {10.1063/1.439715} month = mar, day = "01", } + +@Article{wikipedia_cubic_function, + author = "Wikipedia", + title = "Cubic function", + journal = "Wikipedia", + year = "2009", + month = mar, + day = "23", + url = "http://en.wikipedia.org/wiki/Cubic_equation", +} diff --git a/src/sawsim.nw b/src/sawsim.nw index 03da255..cc50697 100644 --- a/src/sawsim.nw +++ b/src/sawsim.nw @@ -149,8 +149,12 @@ $N$-state rate reactions. This allows simulations of unfolding-refolding, metastable intermediates, etc., but required reorganizing the command line interface. +Version 0.7 added tension model inverses, which seems to reduce +computation time by about a factor of 2. I was expecting more, but +I'll take what I can get. + <>= -#define VERSION "0.6" +#define VERSION "0.7" @ \subsection{License} @@ -354,6 +358,9 @@ typedef struct environment_struct { <>= #ifndef GLOBAL_H #define GLOBAL_H + +#define DOUBLE_PRECISION 1e-12 + <> <> <> @@ -452,6 +459,7 @@ double find_tension(list_t *state_list, environment_t *env, double x, int transi { static int num_active_groups=0; static one_dim_fn_t **per_group_handlers = NULL; + static one_dim_fn_t **per_group_inverse_handlers = NULL; static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */ static double last_x; tension_handler_data_t *tdata; @@ -464,11 +472,11 @@ double find_tension(list_t *state_list, environment_t *env, double x, int transi if (transition != 0 || num_active_groups == 0) { /* setup statics */ /* get new statics, freeing the old ones if they aren't NULL */ - get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_data); + get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data); /* fill in the group handlers and params */ #ifdef DEBUG - fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]); + fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p, (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]); #endif for (i=0; istiffness); + F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness); last_x = x; @@ -704,6 +712,7 @@ model. <>= typedef struct tension_model_getopt_struct { one_dim_fn_t *handler; /* fn implementing the model on a group */ + one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */ char *name; /* name string identifying the model */ char *description; /* a brief description of the model */ int num_params; /* number of extra params for the model */ @@ -712,7 +721,10 @@ typedef struct tension_model_getopt_struct { create_data_func_t *creator; /* param-string -> model-data-struct fn */ destroy_data_func_t *destructor; /* fn to free the model data structure */ } tension_model_getopt_t; -@ +@ Set [[inverse_handler]] to [[NULL]] if you wish to invert using +bisection. Obviously, this will be slower than computing an +analytical inverse and more prone to rounding errors, but it may be +easier if a closed-form inverse does not exist. <>= tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = { @@ -863,8 +875,8 @@ void get_options(int argc, char **argv, #ifdef DEBUG for (i=0; i= 3) tension_models[tension_model].params = strings[num_strings-1]; #ifdef DEBUG - fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s (%s), group %g\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group); + fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s, group %d\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group); #endif INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params); push(pState_list, create_state(state_name, tension_models[tension_model].handler, + tension_models[tension_model].inverse_handler, tension_group, model_params, tension_models[tension_model].destructor, @@ -1311,6 +1324,7 @@ evaluation (see Section \ref{sec.find_tension}). typedef struct state_struct { char *name; /* name string */ one_dim_fn_t *tension_handler; /* tension handler */ + one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */ int tension_group; /* for grouping similar tension states */ void *tension_params; /* pointer to tension parameters */ destroy_data_func_t *destroy_tension; @@ -1336,6 +1350,7 @@ cleaning up. <>= state_t *create_state(char *name, one_dim_fn_t *tension_handler, + one_dim_fn_t *inverse_tension_handler, int tension_group, void *tension_params, destroy_data_func_t *destroy_tension, @@ -1349,6 +1364,7 @@ state_t *create_state(char *name, strcpy(ret->name, name); /* we know ret->name is long enough */ ret->tension_handler = tension_handler; + ret->inverse_tension_handler = inverse_tension_handler; ret->tension_group = tension_group; ret->tension_params = tension_params; ret->destroy_tension = destroy_tension; @@ -1360,7 +1376,7 @@ state_t *create_state(char *name, ret->grouped_states = NULL; #ifdef DEBUG - fprintf(stderr, "generated state %s (%p) with tension handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->tension_params, ret->tension_group); + fprintf(stderr, "generated state %s (%p) with tension handler %p, inverse handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->inverse_tension_handler, ret->tension_params, ret->tension_group); #endif return ret; } @@ -1522,10 +1538,12 @@ active groups. void get_active_groups(list_t *state_list, int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, + one_dim_fn_t ***pPer_group_inverse_handlers, void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */) { void **active_groups=NULL; - one_dim_fn_t *handler, *old_handler, **per_group_handlers=NULL; + one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler; + one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL; void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */ int i, j, num_domains, group, old_group, num_active_groups=0; list_t *active_states_list=NULL; @@ -1541,6 +1559,7 @@ void get_active_groups(list_t *state_list, while (state_list != NULL) { state = S(state_list); handler = state->tension_handler; + inverse_handler = state->inverse_tension_handler; group = state->tension_group; num_domains = state->num_domains; if (list_index(active_states_list, state) == -1) { @@ -1577,6 +1596,8 @@ void get_active_groups(list_t *state_list, /* allocate the arrays we'll be returning */ per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups); assert(per_group_handlers != NULL); + per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups); + assert(per_group_inverse_handlers != NULL); per_group_data = (void **)malloc(sizeof(void *)*num_active_groups); assert(per_group_data != NULL); #ifdef DEBUG @@ -1596,22 +1617,25 @@ void get_active_groups(list_t *state_list, /* populate the arrays we'll be returning */ active_states_list = head(active_states_list); old_handler = NULL; + old_inverse_handler = NULL; j = -1; /* j is the current active _group_ index */ while (active_states_list != NULL) { state = (state_t *) pop(&active_states_list); handler = state->tension_handler; + inverse_handler = state->inverse_tension_handler; group = state->tension_group; num_domains = state->num_domains; - if (handler != old_handler || group != old_group) { + if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) { j++; /* move to the next group */ tdata = (tension_handler_data_t *) per_group_data[j]; per_group_handlers[j] = handler; + per_group_inverse_handlers[j] = inverse_handler; tdata->group_tension_params = NULL; tdata->env = NULL; tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */ } #ifdef DEBUG - fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, group, num_domains); + fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, inverse_handler, group, num_domains); #endif for (i=0; igroup_tension_params, state->tension_params, 1); @@ -1620,12 +1644,15 @@ void get_active_groups(list_t *state_list, fprintf(stderr, "%s:\tpushed %d copies of %p onto data %p's param list %p\n", __FUNCTION__, num_domains, state->tension_params, tdata, tdata->group_tension_params); #endif old_handler = handler; + old_inverse_handler = inverse_handler; old_group = group; } /* free old groups */ if (*pPer_group_handlers != NULL) free(*pPer_group_handlers); + if (*pPer_group_inverse_handlers != NULL) + free(*pPer_group_inverse_handlers); if (*pPer_group_data != NULL) { for (i=0; i<*pNum_active_groups; i++) free((*pPer_group_data)[i]); @@ -1634,10 +1661,11 @@ void get_active_groups(list_t *state_list, /* set new groups */ #ifdef DEBUG - fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]); + fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]); #endif *pNum_active_groups = num_active_groups; *pPer_group_handlers = per_group_handlers; + *pPer_group_inverse_handlers = per_group_inverse_handlers; *pPer_group_data = per_group_data; } @ @@ -2024,7 +2052,8 @@ int main(void) @ <>= -extern double wlc(double x, double T, double p, double L); +<> + START_TEST(test_wlc_at_zero) { double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30; @@ -2172,6 +2201,7 @@ indicates that the groups have changed. <>= double tension_balance(int num_tension_groups, one_dim_fn_t **tension_handlers, + one_dim_fn_t **inverse_tension_handlers, void **params, /* array of pointers to tension_handler_data_t */ double last_x, double x, @@ -2180,12 +2210,13 @@ double tension_balance(int num_tension_groups, <>= double tension_balance(int num_tension_groups, one_dim_fn_t **tension_handlers, + one_dim_fn_t **inverse_tension_handlers, void **params, /* array of pointers to tension_handler_data_t */ double last_x, double x, double *stiffness) { - static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL}; + static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL}; double F, xo_guess, xo, lb, ub; double min_relx=1e-6, min_rely=1e-6; int max_steps=200, i; @@ -2202,9 +2233,10 @@ double tension_balance(int num_tension_groups, /* construct new x_xo_data */ x_xo_data.n_groups = num_tension_groups; x_xo_data.tension_handler = tension_handlers; + x_xo_data.inverse_tension_handler = inverse_tension_handlers; x_xo_data.handler_data = params; #ifdef DEBUG - fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0]); + fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p), x_of_xo_data %p\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0], &x_xo_data); #endif x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups); assert(x_xo_data.xi != NULL); @@ -2280,6 +2312,7 @@ double tension_balance(int num_tension_groups, typedef struct x_of_xo_data_struct { int n_groups; one_dim_fn_t **tension_handler; /* array of fn pointers */ + one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */ void **handler_data; /* array of pointers to tension_handler_data_t */ double *xi; /* array of group extensions */ } x_of_xo_data_t; @@ -2296,7 +2329,7 @@ double x_of_xo(double xo, void *pdata) assert(data->n_groups > 0); data->xi[0] = xo; #ifdef DEBUG - fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0]); + fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p (x_xo_data %p)\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0], data); fflush(stderr); #endif F = (*data->tension_handler[0])(xo, data->handler_data[0]); @@ -2307,19 +2340,22 @@ double x_of_xo(double xo, void *pdata) 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]; + if (data->inverse_tension_handler[i] != NULL) { + xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]); + } else { /* invert numerically */ + if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */ + else xi_guess = data->xi[i]; #ifdef DEBUG - fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]); + fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]); #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_relx, min_rely, max_steps, NULL); + 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_relx, min_rely, max_steps, NULL); + } + x += xi; + data->xi[i] = xi; #ifdef DEBUG - fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F); + fprintf(stderr, "%s:\tfound 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:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x); @@ -2444,6 +2480,9 @@ void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double { double yx, step, x=xguess; int i=0; +#ifdef DEBUG + fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params); +#endif yx = (*f)(x, params); #ifdef DEBUG fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params); @@ -2542,8 +2581,9 @@ START_TEST(test_single_function) { double k=5, x=3, last_x=2.0, F; one_dim_fn_t *handlers[] = {&hooke}; + one_dim_fn_t *inverse_handlers[] = {NULL}; void *data[] = {&k}; - F = tension_balance(1, handlers, data, last_x, x, NULL); + F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL); fail_unless(F = k*x, NULL); } END_TEST @@ -2556,8 +2596,9 @@ START_TEST(test_double_hooke) { double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e; one_dim_fn_t *handlers[] = {&hooke, &hooke}; + one_dim_fn_t *inverse_handlers[] = {NULL, NULL}; void *data[] = {&k1, &k2}; - F = tension_balance(2, handlers, data, last_x, x, NULL); + F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL); x1e = x*k2/(k1+k2); Fe = k1*x1e; x2e = x1e*k1/k2; @@ -3557,12 +3598,17 @@ details For unstretchable domains. <>= -{NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL} +{NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL} @ \subsection{Constant} \label{sec.const_tension} +An infinitely stretchable domain providing a constant tension. +Obviously this cannot be inverted, so you can't put this domain in +series with any other domains. We include it mostly for testing +purposes. + <>= <> <> @@ -3575,8 +3621,6 @@ For unstretchable domains. @ -An infinitely stretchable domain providing a constant tension. - <>= extern double const_tension_handler(double x, void *pdata); @ @@ -3651,14 +3695,16 @@ char const_tension_param_string[] = "0"; @ <>= -{&const_tension_handler, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t} +{&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t} @ \subsection{Hooke} \label{sec.hooke} <>= -<> +<> +<> +<> <> @ @@ -3678,16 +3724,17 @@ For a simple proof, see Appendix \ref{sec.math_hooke}. <>= extern double hooke_handler(double x, void *pdata); +extern double inverse_hooke_handler(double F, void *pdata); @ -<>= -double hooke_handler(double x, void *pdata) -{ - tension_handler_data_t *data = (tension_handler_data_t *)pdata; +First we define a function that computes the effective spring constant +(stored in a single [[hooke_param_t]]) for the entire group. +<>= +static void hooke_param_grouper(tension_handler_data_t *data, + hooke_param_t *param) { list_t *list = data->group_tension_params; double k=0.0; - assert (x >= 0.0); list = head(list); assert(list != NULL); /* empty active group?! */ while (list != NULL) { @@ -3704,19 +3751,51 @@ double hooke_handler(double x, void *pdata) fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n", __FUNCTION__, k, x, k*x, data); #endif - return k*x; + param->k = k; +} + +@ + +Everyone knows Hooke's law: $F=kx$. +<>= +double hooke_handler(double x, void *pdata) +{ + hooke_param_t param = {0}; + + assert (x >= 0.0); + hooke_param_grouper((tension_handler_data_t *)pdata, ¶m); + return param.k*x; +} + +@ + +Inverting Hooke's law gives $x=F/k$. +<>= +double inverse_hooke_handler(double F, void *pdata) +{ + hooke_param_t param = {0}; + + assert (F >= 0.0); + hooke_param_grouper((tension_handler_data_t *)pdata, ¶m); + return F/param.k; } + @ +Not much to keep track of for a single effective spring, just the +spring constant $k$. <>= typedef struct hooke_param_struct { double k; /* spring constant in N/m */ } hooke_param_t; + @ +We wrap up our Hooke implementation with some book-keeping functions. <>= extern void *string_create_hooke_param_t(char **param_strings); extern void destroy_hooke_param_t(void *p); + @ <>= @@ -3753,7 +3832,7 @@ char hooke_param_string[]="0.05"; @ <>= -{hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t} +{hooke_handler, inverse_hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t} @ \subsection{Worm-like chain} @@ -3761,18 +3840,27 @@ char hooke_param_string[]="0.05"; We can model several unfolded domains with a single worm-like chain. <>= -<> +<> <> +<> <> @ <>= + <> +<> <> <> @ +<>= +<> +<> +<> +@ + The combination of all unfolded domains is modeled as a worm like chain, with the force $F_{\text{WLC}}$ approximately given by $$ @@ -3785,7 +3873,7 @@ $T$ is the absolute temperature, $p$ is the persistence length, and $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}. <>= -double wlc(double x, double T, double p, double L) +static double wlc(double x, double T, double p, double L) {/* N m K m m */ double xL=0.0; /* uses global k_B */ assert(x >= 0); @@ -3796,19 +3884,82 @@ double wlc(double x, double T, double p, double L) // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL )); return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ); } + @ -This model assumes that each unfolded domain has the same persistence length. -<>= -extern double wlc_handler(double x, void *pdata); +Inverting the approximate WLC is fairly straightforward, if a bit tedious. +Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have +\begin{align} + F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\ + F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\ + 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\ + &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2}) + + x_L - 2x_L^2 + x_L^3 \\ + &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}] + + x_L \p[{2F_T + \frac{1}{2} + 1}] + + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}] + + x_L^3 \\ + &\approx -F_T + + x_L \p({2F_T + \frac{3}{2}}) + - x_L^2 \p({F_T + \frac{9}{4}}) + + x_L^3 \\ + &\approx ax_L^3 + bx_L^2 + cx_L + d +\end{align} +where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and +$d \equiv -F_T$. +% From \citet{wikipedia_cubic_function} the discriminant +%\begin{align} +% \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\ +% &= -4F_T\p({F_T + \frac{9}{4}})^3 +% + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2 +% - 4 \p({2F_T + \frac{3}{2}})^3 +% + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}}) +% - 27F_T^2 \\ +% &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}}) +%% a^3 + 3a^2b + 3ab^2 + b^3 +% + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}}) +% \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\ +% &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}}) +%% a^3 + 3a^2b + 3ab^2 + b^3 +% + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}}) +% - 27F_T^2 \\ +% &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T +% + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\ +% &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2} +% + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T +% - 27F_T^2 \\ +% &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36}) +% + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\ +% &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}}) +% \p({\frac{729}{64} - \frac{27}{2}}) \\ +% &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64} +%\end{align} +We can plug these coefficients into the GSL cubic solver to invert the +WLC\citep{galassi05}. The applicable root is always the one which +comes before the singularity, which will be the smallest real root. +<>= +static double inverse_wlc(double F, double T, double p, double L) +{/* m N K m m */ + double FT = F*p/(k_B*T); /* uses global k_B */ + double xL0, xL1, xL2; + int num_roots; + assert(F >= 0); + assert(T > 0); assert(p > 0); assert(L > 0); + num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2); + assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1); + if (xL0 < 0) xL0 = 0.0; + return xL0*L; +} + @ -<>= -double wlc_handler(double x, void *pdata) -{ - tension_handler_data_t *data = (tension_handler_data_t *)pdata; +First we define a function that computes the effective WLC parameters +(stored in a single [[wlc_param_t]]) for the entire group. +<>= +static void wlc_param_grouper(tension_handler_data_t *data, + wlc_param_t *param) { list_t *list = data->group_tension_params; - double p, L=0.0; + double p=0.0, L=0.0; list = head(list); assert(list != NULL); /* empty active group?! */ @@ -3827,8 +3978,44 @@ double wlc_handler(double x, void *pdata) 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); + param->p = p; + param->L = L; } + +@ + +<>= +extern double wlc_handler(double x, void *pdata); +@ + +This model requires that each unfolded domain in the group have the +same persistence length. +<>= +double wlc_handler(double x, void *pdata) +{ + tension_handler_data_t *data = (tension_handler_data_t *)pdata; + wlc_param_t param = {0}; + + wlc_param_grouper(data, ¶m); + return wlc(x, data->env->T, param.p, param.L); +} + +@ + +<>= +extern double inverse_wlc_handler(double F, void *pdata); +@ + +<>= +double inverse_wlc_handler(double F, void *pdata) +{ + tension_handler_data_t *data = (tension_handler_data_t *)pdata; + wlc_param_t param = {0}; + + wlc_param_grouper(data, ¶m); + return inverse_wlc(F, data->env->T, param.p, param.L); +} + @ <>= @@ -3863,6 +4050,7 @@ void destroy_wlc_param_t(void *p /* wlc_param_t * */) if (p) free(p); } + @ We define [[wlc_param_string]] as a global array because hard-coding the values into the tension model getopt structure gives a read-only version. @@ -3889,7 +4077,7 @@ char wlc_param_string[]="0.39e-9,28.4e-9"; @ <>= -{wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t} +{wlc_handler, inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t} @ Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696. @@ -4019,7 +4207,7 @@ char fjc_param_string[]="0.5e-9,200"; @ <>= -{fjc_handler, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t} +{fjc_handler, NULL, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t} @ \subsection{Tension model implementation} @@ -5628,8 +5816,8 @@ the remaining two [[pdflatex]] passes insert the citations and finalize the vari semi-clean_latex : rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \ $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \ - $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \ - $(BUILD_DIR)/sawsim.bib + $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \ + $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib clean_latex : semi-clean_latex rm -f $(DOC_DIR)/sawsim.pdf @ diff --git a/testing/bell_rate/run_test.sh b/testing/bell_rate/run_test.sh index 931eeba..78421ca 100755 --- a/testing/bell_rate/run_test.sh +++ b/testing/bell_rate/run_test.sh @@ -26,13 +26,17 @@ # The histogram scaling isn't dD/df, though, it's dD/d(bin) # dD/d(bin) = dD/df df/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})} # -# test with lots of runs of v-dx=1 -# sawsim -v1 -mhooke -a1 -kbell -K1,1 -T7.24296e22 -F1 -# ^- v=1 ^- k=1 ^-a=1 ^-1/kB +# test with lots of runs of +# v- T=1/kB v- v=1 v- k=1 +# sawsim -T7.24296e22 -v1 -s folded,hooke,1 -N1 -s unfolded,null \ +# -k "folded,unfolded,bell,{1,1}" -q folded +# ^- {a,dx}={1,1} # (kB = 1.3806503 10-23 J/K, we use T=1/kB to get the bell K=e^{-f}) # -# 70 runs of that takes ~1 second: -# time xtimes 200 sawsim -v1 -mhooke -a1 -kconst -K1 -F1 > /dev/null +# 60 runs of that takes ~1 second: +# time xtimes 60 sawsim -T7.24296e22 -v1 -s folded,hooke,1 -N1 \ +# -s unfolded,null -k "folded,unfolded,bell,{1,1}" -q folded \ +# > /dev/null # # usage: bell_rate.sh @@ -49,6 +53,7 @@ PNGFILE=bell_rate.png # dD/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})} # where z ≡ a/kv, and b = Δx/kBT KB=1.3806503e-23 # J/K +#STYLE="ones" STYLE="protein" if [ "$STYLE" == "ones" ] then @@ -74,9 +79,15 @@ THEORY_FN="$df * $N * $Z * exp($B*x + $Z/$B *(1 - exp($B*x)))" > $DATA if [ "$ISACLUSTER" -eq 1 ] then - qcmd "xtimes $N $SAWSIM -v$V -mhooke -a$K -kbell -K$A,$DX -T$T -F1 | grep -v '^#' >> $DATA" + # Sawsim <= 0.5 + #qcmd "xtimes $N $SAWSIM -T$T -v$V -mhooke -a$K -kbell -K$A,$DX -F1 | grep -v '^#' >> $DATA" + # Sawsim >= 0.6 + qcmd "xtimes $N $SAWSIM -T$T -v$V -s folded,hooke,$K -N1 -s unfolded,null -k 'folded,unfolded,bell,{$A,$DX}' -q folded | grep -v '^#' | cut -f1 >> $DATA" else - xtimes $N $SAWSIM -v$V -mhooke -a$K -kbell -K$A,$DX -T$T -F1 | grep -v '^#' >> $DATA + # Sawsim <= 0.5 + #xtimes $N $SAWSIM -T$T -v$V -mhooke -a$K -kbell -K$A,$DX -F1 | grep -v '^#' >> $DATA + # Sawsim >= 0.6 + xtimes $N $SAWSIM -T$T -v$V -s folded,hooke,$K -N1 -s unfolded,null -k "folded,unfolded,bell,{$A,$DX}" -q folded | grep -v '^#' | cut -f1 >> $DATA fi # histogram the data @@ -87,7 +98,7 @@ FIT=`cat $HIST | ./fit_bell` #echo "$FIT" FITZ=`echo "$FIT" | sed -n 's/a:[[:space:]]*//p'` FITB=`echo "$FIT" | sed -n 's/b:[[:space:]]*//p'` -echo -e "a:\t$FITZ\t(expected $Z)" +echo -e "z:\t$FITZ\t(expected $Z)" echo -e "b:\t$FITB\t(expected $B)" # generate gnuplot script diff --git a/testing/common/fit.pyc b/testing/common/fit.pyc deleted file mode 100644 index 396d9d25504bf2b1931414ce0d9ad0fff0a2b2ea..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 4361 zcma)9ZFAek5k8Q5F%`LS?9^j7NqaMCVydNRDoz_$88?k3*Vcz94%Dd|r5+9h97$Ms zp>YR^^4R?l=VO0Ue^mPi`aFA($hec4Qo`{7-0kf?yZh`i|M%s!-T(aSQLLIz8{dDz zWe+eUN*!Vh)N!B=TQX^#Pul9RtrFZWsidvUvKj>{SyDeJ^;)TK4_8#O%$vhim8{6* z689f;S1J4@hCioHPg3j3aTVn{iCm=PGEa-KP77TY+PVZ2SL#UD)>NB`8Ks3ubY3PV z>uddOS+ChiS!anptx{*OZhGb>+!`6EPeo3Fki5s`5p(D z?{Vm;WK{rJQ{h4P5+x7YfTOH5tWWdwMWe7bRaI8ox-zb=3aj13=%TFhC`(@ehbp%^ zNk=17nIbl88lNNMPK_ym1mU#C>8MC>H>!)+r9kGa`?n21w+Dc$#8gBV*g|q@ML?VZ zubw_DYn|724xxsxD5^X$l?U77D8A^Sa9vgdSap*6gjn<9a#ZW1{W8fCVe z3z5}c(6~0%4q?7aqipEKbvJr)DC{m^xYJW8hv1Mkf|3aEi)qTDPBil`FiTHNHnrLm zwyumGmnHD1U71IM$Ed7CMvd1~Q=O!i+8XQN0tnbOCnhVW^oQ)6n7Tp`rLooLI~Ua7 z^O`*tEV>&6p|*9DWoP2LNd(A&DfD#`u@oE0@vyLjm(w`P7E*$G^w(qogX5?yYOAMZ zTDZnVaRfiYBUPH1#`tN4*P;^?fX8IF*HK=bE7v1lWb3vCe<@LOda zeGcJX%W|Yo;mBJSU{Lq6U?}k(ap!zA9T^RUj!Fa-_r{rtDjyZ0$xVJ_Dx6HW;Ggl_ zV6Bg;sE8*bPZ!1{cw0iOwq=$kczh6@_z)L$PfMyds|n>);ScY}ME z)l0t?Dc=UaEm8IuLlX|jRnbYDfaoQtvC*fgn;~UDzz`%>Zkza$;6s$1jC?Goy+$NFe7*DJ2KBv zTAZedo45{ciaKw?gBwEfUP!hMxxS0b5AF(P{r+vHcUjy9v-ZKels1$E9u9wrxn0Ac zS}N$YZni!QVoF=VMSd0dBTP^bUbWTib3wVKcApN^-&^X{vU&yJJ@$ZC5PscKH(tS8 z;v|57?KA9renmaO6C4f$=$n>;NX?be0EaYwxIKz@Iaqt*lLoppcALOEHMda5a41WZu!y|@F%m&&nLz(zsg+H)wC`1>CSuE*7BFb$HRChf zS@H)|_a@LnJ`fYWz+`rpJrH8oGp0{B#3ekx`}93YUnI3ON94tp8nx8*hmW4Wub%vW zh8JmUP%#VSQo=MtrEaqABEgzqv5b9rIQ%H--0OTaUfcJXL5;=fs48OW{ebUDT5qDm09h!T+qkutFc6;ikvSJdg~2Eg{RR`@1G|BY&le=mbXAh&d*%#U z#6W!yfEWf;)biwc80?Yu(KkuMgqD!hu;5}a%mh6Q+o8tAqxM#R&=`JUsSjMTx1-@- zTErAPo(F6xg_h>puuO^BBLdo=3sDu{S+X?>2$X zayYl}aaEqMWAPzgIHX=Wrok8O6^s8gjgOo5h~$g12?G$^(f~dL5W>r&{Q+V;Be@eBf}x2*BNgK>TsBp@o*S2sf9MecRBk7hxa)A ziUU!T6XtNrOcK z4v4cIAxR_7o?jMjM)?F+W5W{28jC>Gues@(5H!*kd@VtG_4$4E4=I}5`ULx$urxR@ zG*Pn9k+K}>QB~$T0w>`&mz_qj*|Y||bL30{?mjIy<1(xBLYVu@gPpxSG{SIztmpUUM~K-KIFmr{;(I@_R@c4#An5zCL;;q@3Z`pVY``I>pZlc{coR>#1k= zX_h)Ed?AZlm(YbDbJ*Y@fqjcJM#9`3zBXRMvkI2(!n{YAFVx3XwGzCeR)f3q4yKJ2 zJ6tuwkSKKwr=WWjvhxiHyBA3TfXRv|JxsbPl@lIA$?R|*aJEM0fGYxck@s(rE#Xbh z*JZAnaXEyvTZ@5l=J$oP-zmgj}{B-h`6?Rdw+VzC4oPinsB%9;^kIF-rXoe~s6| zO;VWq=S4ZZj^}cURQTij;_~w78TtcB88I>3!YZ%Y9j@oCe` Qt+cM9JzZ&Czt`#f7iD_5RsaA1 diff --git a/testing/const_rate/const_rate.sh b/testing/const_rate/const_rate.sh index 5b59829..56ae14f 100755 --- a/testing/const_rate/const_rate.sh +++ b/testing/const_rate/const_rate.sh @@ -54,15 +54,19 @@ THEORY_FN="$theoryA * exp(x / ($theoryB))" > $DATA if [ "$ISACLUSTER" -eq 1 ] then - if [ $DEBUG -eq 1 ]; then - echo "qcmd \"xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA\"" - fi - qcmd "xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA" + # Sawsim <= 0.5 + #qcmd "xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA" + # Sawsim >= 0.6 + qcmd xtimes $N $SAWSIM -v$V -s folded,hooke,$K -N1 -s unfolded,null \ + -k folded,unfolded,const,$RATE -q folded \ + | grep -v '^#' | cut -f1 >> $DATA else - if [ $DEBUG -eq 1 ]; then - echo "xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA" - fi - xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA + # Sawsim <= 0.5 + #xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA + # Sawsim >= 0.6 + xtimes $N $SAWSIM -v$V -s folded,hooke,$K -N1 -s unfolded,null \ + -k folded,unfolded,const,$RATE -q folded \ + | grep -v '^#' | cut -f1 >> $DATA fi # histogram the data -- 2.26.2