From: W. Trevor King Date: Wed, 25 Feb 2009 22:27:23 +0000 (-0500) Subject: Added stiffness env. parameter and stiffness-corrected Bell model. X-Git-Tag: v0.5 X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=6b92bcd6412cbd40ae1c76ea23f6daeab332939c;p=sawsim.git Added stiffness env. parameter and stiffness-corrected Bell model. Following Walton et al. and my poster. Currently a fairly ugly hack to get some simulation data for my poster (in 4 days!). I will come back and clean things up afterwards. Due to the coarseness of the tension balancer, it's probably a better idea to get stiffnesses for each of the groups separately and then add the stiffnesses in series. This removes the balancer completely. --- diff --git a/src/sawsim.nw b/src/sawsim.nw index 9a41951..1253c54 100644 --- a/src/sawsim.nw +++ b/src/sawsim.nw @@ -1,5 +1,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%% -% we have our own preamble, +% we have our own preamble, % use `noweave -delay` to not print noweb's automatic one % -*- mode: Noweb; noweb-code-mode: c-mode -*- \documentclass[letterpaper, 11pt]{article} @@ -133,9 +133,12 @@ be controlled from the command line. Also added interpolating lookup tables to accelerate slow unfolding rate model evaluation, and command line control of tension groups. +Version 0.5 added the stiffness environmental parameter and it's +associated unfolding models. + <>= -#define VERSION "0.4" -@ +#define VERSION "0.5" +@ \subsection{License} @@ -150,7 +153,7 @@ line control of tension groups. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License @@ -161,14 +164,14 @@ line control of tension groups. The author may be contacted at on the Internet, or write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St., Philadelphia PA 19104, USA. -@ +@ <>= /* <> */ -@ +@ \section{Main} @@ -195,6 +198,7 @@ function pointers, with implementations for \begin{packed_item} \item Null (Appendix \ref{sec.null_k}), \item Bell's model (Appendix \ref{sec.bell}), + \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}), \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and \item Kramers' saddle point model (Appendix \ref{sec.kramers}). \end{packed_item} @@ -219,7 +223,7 @@ int main(int argc, char **argv) get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded); setup(); - if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n"); + if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n"); if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n"); while (num_folded > 0) { F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding); @@ -228,7 +232,7 @@ int main(int argc, char **argv) else dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, 0); unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded); - if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F); + if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness); if (xstep == 0) { x += v*dt; } else { @@ -267,11 +271,12 @@ calculated, bridging the gap between theory and experiment. Environmental parameters important in determining reaction rates and tensions (e.g. temperature, pH) are stored in a single structure to facilitate extension to more complicated models in the future. -<>= +<>= typedef struct environment_struct { double T; + double stiffness; } environment_t; -@ +@ <>= #ifndef GLOBAL_H @@ -281,7 +286,7 @@ typedef struct environment_struct { <> <> #endif /* GLOBAL_H */ -@ +@ Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being included multiple times. @@ -293,7 +298,7 @@ included multiple times. <> <> <> -@ +@ \subsection{Tension} \label{sec.find_tension} @@ -334,13 +339,13 @@ grouping domains with seperate models doesn't make sense. <>= #define NUM_TENSION_MODELS 5 -@ +@ <>= extern one_dim_fn_t *tension_handlers[]; @ -<>= +<>= one_dim_fn_t *tension_handlers[] = { NULL, &const_tension_handler, @@ -349,11 +354,11 @@ one_dim_fn_t *tension_handlers[] = { &fjc_handler, }; -@ +@ <>= <> -@ +@ <>= typedef struct tension_handler_data_struct { @@ -364,7 +369,7 @@ typedef struct tension_handler_data_struct { void *persist; } tension_handler_data_t; -@ +@ After sorting the chain into separate groups $G_i$, with tension handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of @@ -373,7 +378,11 @@ each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\; class. [[find_tension]] is basically a wrapper to feed proper input into the [[tension_balance]] function. [[unfolding]] should be set to 0 if there were no unfoldings since the previous call to -[[find_tension]]. +[[find_tension]]. While were messing with the tension handlers, we +also grab a measure of the current linker stiffness for the +environmental variable, which is needed by some of the unfolding rate +models (e.g. the linker-stiffness corrected Bell model, Appendix +\ref{sec.kbell}). <>= double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding) { @@ -390,7 +399,7 @@ double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, l list_print(stderr, domain_list, "domain list"); #endif - if (unfolding != 0 || num_active_groups == 0) { /* setup statics */ + if (unfolding != 0 || num_active_groups == 0) { /* setup statics */ /* free old statics */ if (per_group_handlers != NULL) free(per_group_handlers); @@ -419,7 +428,32 @@ double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, l } /* else roll with the current statics */ F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x); + + { + double dx; + double xlow; + double Flow; +#ifdef DEBUG + fprintf(stderr, "finding linker stiffness at distance %g, tension %g\n", x, F); +#endif + if (last_x == -1.0) + dx = x/200.0; + else + dx = (x-last_x); + if (dx == 0) { /* e.g. if x == 0 */ + env->stiffness = 0; + } else { + xlow = x-dx; + Flow = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, xlow); + env->stiffness = (F-Flow)/dx; + } +#ifdef DEBUG + fprintf(stderr, " stiffness (%g-%g)/%g = %g\t(dx = %g = %g-%g)\n", F, Flow, dx, env->stiffness, dx, x, xlow); +#endif + } + last_x = x; + return F; } @ For the moment, we will restrict our group boundaries to a single @@ -490,12 +524,12 @@ adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim shouldn't be a problem. To reassure yourself, you can ask the simulator to print the smallest timestep that was used with TODO. <>= -int random_unfoldings(list_t *domain_list, double tension, +int random_unfoldings(list_t *domain_list, double tension, double dt, environment_t *env, int *num_folded_domains) { /* returns 1 if there was an unfolding and 0 otherwise */ while (domain_list != NULL) { - if (D(domain_list)->state == DS_FOLDED + if (D(domain_list)->state == DS_FOLDED && domain_unfolds(tension, dt, env, D(domain_list))) { if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) fprintf(stdout, "%g\n", tension); @@ -507,7 +541,7 @@ int random_unfoldings(list_t *domain_list, double tension, } return 0; } -@ +@ \subsection{Adaptive timesteps} \label{sec.adaptive_dt} @@ -528,7 +562,7 @@ below a set threshold, so I've removed it to Appendix \ref{app.adaptive_dt}. \section{Models} - + TODO: model intro. The models provide the physics of the simulation, but the simulation @@ -541,7 +575,7 @@ and unfolding models are defined in Appendix \ref{sec.k_models}. <>= #define k_B 1.3806503e-23 /* J/K */ -@ +@ \section{Command line interface} @@ -552,7 +586,7 @@ and unfolding models are defined in Appendix \ref{sec.k_models}. <> <> <> -@ +@ \subsection{Model selection} \label{app.model_selection} @@ -564,17 +598,17 @@ of available models, containing their important parameters. <>= <> -@ +@ <>= typedef void *create_data_func_t(char **param_strings); typedef void destroy_data_func_t(void *param_struct); -@ +@ <>= <> <> -@ +@ \subsubsection{Tension} @@ -593,7 +627,7 @@ 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; -@ +@ <>= tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = { @@ -603,12 +637,12 @@ tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = { <>, <> }; -@ +@ \subsubsection{Unfolding rate} <>= -#define NUM_K_MODELS 5 +#define NUM_K_MODELS 6 typedef struct k_model_getopt_struct { char *name; @@ -620,17 +654,18 @@ typedef struct k_model_getopt_struct { create_data_func_t *creator; destroy_data_func_t *destructor; } k_model_getopt_t; -@ +@ <>= k_model_getopt_t k_models[NUM_K_MODELS] = { <>, <>, <>, + <>, <>, <> }; -@ +@ \subsection{help} @@ -683,7 +718,7 @@ void help(char *prog_name, double P, double dt_max, double v, double xstep, printf(" 123.456e-12\n"); printf(" ...\n"); printf("In verbose mode, the entire Force-vs-distance curve is output:\n"); - printf(" #Position (m)\tForce (N)\n"); + printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n"); printf(" 0.001\t0.0023\n"); printf(" ...\n"); printf("-V\tChange output to verbose mode\n"); @@ -709,7 +744,7 @@ void help(char *prog_name, double P, double dt_max, double v, double xstep, printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n"); printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name); } -@ +@ \subsection{Get options} @@ -840,7 +875,7 @@ void get_options(int argc, char **argv, assert(env->T > 0.0); return; } -@ +@ <>= int index_tension_model(int n_models, tension_model_getopt_t *models, char *name) @@ -855,7 +890,7 @@ int index_tension_model(int n_models, tension_model_getopt_t *models, char *name } return i; } -@ +@ <>= int index_k_model(int n_models, k_model_getopt_t *models, char *name) @@ -870,7 +905,7 @@ int index_k_model(int n_models, k_model_getopt_t *models, char *name) } return i; } -@ +@ <>= /* requires int num_param_args and char **param_args in the current scope @@ -897,7 +932,7 @@ int index_k_model(int n_models, k_model_getopt_t *models, char *name) param_pointer = NULL; \ free_parsed_list(num_param_args, param_args); \ } while (0); -@ +@ <>= void *generate_domain(enum domain_state_t state, @@ -929,7 +964,7 @@ void *generate_domain(enum domain_state_t state, unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor ); } -@ +@ \phantomsection \appendix @@ -947,7 +982,7 @@ The general layout of our simulation code is: <> <> <
> -@ +@ We include [[math.h]], so don't forget to link to the libm with `-lm'. <>= @@ -965,7 +1000,7 @@ We include [[math.h]], so don't forget to link to the libm with `-lm'. #include "tension_model.h" #include "parse.h" #include "accel_k.h" -@ +@ <>= <> @@ -975,12 +1010,12 @@ We include [[math.h]], so don't forget to link to the libm with `-lm'. <> <> <> -@ +@ <>= <> <> -@ +@ <>= <> @@ -988,19 +1023,19 @@ We include [[math.h]], so don't forget to link to the libm with `-lm'. <> <> <> -@ +@ <>= <> <> -@ +@ <>= void setup(void) { srand(getpid()*time(NULL)); /* seed rand() */ } -@ +@ <>= /* in octal b/c of prefixed '0' */ @@ -1018,7 +1053,7 @@ static unsigned long int flags = 0; <>= #define MAX(a,b) ((a)>(b) ? (a) : (b)) #define MIN(a,b) ((a)<(b) ? (a) : (b)) -@ +@ Note that [[STRMATCH]] chokes if one of the strings is [[NULL]]. <>= @@ -1029,7 +1064,7 @@ static char *temp_string_B; strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \ !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) ) /* +1 to also compare the '\0' */ -@ +@ We also define a macro for our [[check]] unit testing <>= @@ -1044,7 +1079,7 @@ We also define a macro for our [[check]] unit testing -(received-expected)/expected, max_err, #received, \ expected, received); \ } while(0) -@ +@ <>= int happens(double probability) @@ -1069,7 +1104,7 @@ In the less common, domains did unfold in the last step, and we expect the next step to be considerably longer than the previous one. <>= double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, - list_t *domain_list, + list_t *domain_list, environment_t *env, double x, double target_prob, double max_dt, double v) { /* Returns the timestep dt in seconds for the current folded domain. @@ -1087,7 +1122,7 @@ double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, //printf("x %g\tF %g\tk %g\n", x, F, k); dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */ if (dtU > max_dt) { - //printf("overshot max_dt\n"); + //printf("overshot max_dt\n"); dtU = max_dt; } if (v == 0) /* discrete stepping, so force is temporatily constant */ @@ -1108,14 +1143,14 @@ double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, //printf("Try for dt = %g (F = %g)\n", dt, F); k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params); /* returned k may be -1.0 */ - //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); + //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); while (k == -1.0) { /* reduce step until we hit a valid k */ dtU /= 2.0; dt = dtU; /* hopefully, we can use the max dt, see if it works */ F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0); //printf("Try for dt = %g (F = %g)\n", dt, F); k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params); - //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); + //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); } assert(dtU > 1e-14); /* timestep to valid k too small */ dtUCur = target_prob / k; /* safe timestep back from x+dtU */ @@ -1133,14 +1168,14 @@ double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, if (dtCur > dt) /* safe timestep back from x+dt covers dt */ dtL = dt; else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */ - dtU = dt; /* ... stepping out only dtCur would be safe */ + dtU = dt; /* ... stepping out only dtCur would be safe */ dtUCur = dtCur; } else break; /* dtCur = dt */ } return MAX(dtUCur, dtL); } -@ +@ To determine $dt$ for an array of potentially different folded domains, we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains. @@ -1155,7 +1190,7 @@ double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, l double dt=dt_max, new_dt; assert(target_prob > 0.0); assert(target_prob < 1.0); assert(dt_max > 0.0); - + while (domain_list != NULL) { if (D(domain_list)->state == DS_FOLDED) { new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v); @@ -1165,7 +1200,7 @@ double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, l } return dt; } -@ +@ \subsection{Domain data} @@ -1202,7 +1237,7 @@ typedef struct domain_struct { #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \ ? ((domain_t *)(list)->d)->folded_params \ : ((domain_t *)(list)->d)->unfolded_params) -@ +@ [[k]] is a pointer to the function determining the unfolding rate for a given tension. [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]]. [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension. @@ -1283,7 +1318,7 @@ void destroy_domain_list(list_t *domain_list) <> <> <> -@ +@ <>= one_dim_fn_t *get_tension_handler(domain_t *domain) @@ -1297,7 +1332,7 @@ one_dim_fn_t *get_tension_handler(domain_t *domain) return domain->unfolded_handler; } } -@ +@ <>= int get_group(domain_t *domain) @@ -1311,7 +1346,7 @@ int get_group(domain_t *domain) return domain->unfolded_group; } } -@ +@ We already know all possible tension classes, since they are all hardcoded into \prog. However, there may be any number of possible @@ -1349,7 +1384,7 @@ list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) { #endif return ret; } -@ +@ Search a [[list]] of domains, and determine whether a particular model class and group number combination is active. @@ -1364,7 +1399,7 @@ int is_group_active(list_t *list, one_dim_fn_t *handler, int group) } return 0; } -@ +@ Search a [[list]] of domains, and return all domains matching a particular model class and group number. @@ -1384,7 +1419,7 @@ list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group) } return ret; } -@ +@ Because all the node data in lists returned by [[get_group_list]] is also in the main domain list, you shouldn't destroy the node data popped off when destroying the group lists. It will all get cleaned @@ -1421,7 +1456,7 @@ void get_active_groups(list_t *list, } } - /* allocate the array we'll be returning */ + /* allocate the array we'll be returning */ num_active_groups = list_length(active_groups_list); active_groups = (void **) malloc(sizeof(void *)*num_active_groups); assert(active_groups != NULL); @@ -1448,7 +1483,7 @@ void get_active_groups(list_t *list, *pPer_group_handlers = per_group_handlers; *pActive_groups = active_groups; } -@ +@ \section{String parsing} @@ -1466,23 +1501,23 @@ We implement this parsing in [[parse.c]], define the interface in [[parse.h]], a <> <> #endif /* PARSE */ -@ +@ <>= NW_SAWSIM_MODS += parse CHECK_BINS += check_parse -check_parse_MODS = -@ +check_parse_MODS = +@ <>= #define SEP ',' /* argument separator character */ -@ +@ <>= extern void parse_list_string(char *string, char sep, char deeper, char shallower, int *num, char ***string_array); extern void free_parsed_list(int num, char **string_array); -@ +@ [[parse_list_string]] allocates memory, don't forget to free it afterward with [[free_parsed_list]]. It does not alter the original. @@ -1572,7 +1607,7 @@ void free_parsed_list(int num, char **string_array) <> #include "parse.h" <> -@ +@ <>= #include /* assert() */ @@ -1580,7 +1615,7 @@ void free_parsed_list(int num, char **string_array) #include /* fprintf(), stdout *//*!!*/ #include /* strlen() */ #include "parse.h" -@ +@ \subsection{Parsing unit tests} @@ -1591,7 +1626,7 @@ Here we check to make sure the various functions work as expected, using \citeta <> <> <
> -@ +@ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE, atof() */ @@ -1600,12 +1635,12 @@ Here we check to make sure the various functions work as expected, using \citeta #include /* strlen() */ <> #include "parse.h" -@ +@ <>= <> <> -@ +@ <>= Suite *test_suite (void) @@ -1616,7 +1651,7 @@ Suite *test_suite (void) <> return s; } -@ +@ <>= @@ -1699,14 +1734,14 @@ Here we check to make sure the various functions work as expected, using \citeta <> <> <
> -@ +@ <>= #include -@ +@ <>= -@ +@ <>= <> @@ -1715,7 +1750,7 @@ Here we check to make sure the various functions work as expected, using \citeta <> <> <> -@ +@ <>= Suite *test_suite (void) @@ -1740,7 +1775,7 @@ Suite *test_suite (void) */ return s; } -@ +@ <
>= int main(void) @@ -1753,18 +1788,18 @@ int main(void) srunner_free(sr); return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE; } -@ +@ \subsection{F tests} <>= <> -@ +@ <>= <> -@ +@ <>= <> -@ +@ <>= extern double wlc(double x, double T, double p, double L); @@ -1785,23 +1820,23 @@ START_TEST(test_wlc_at_half) * = 0.25 * 3 = 3/4 * linear term = x/L = 0.5 * nonlinear + linear = 0.75 + 0.5 = 1.25 - * wlc = 10e21*1.25 = 12.5e21 + * wlc = 10e21*1.25 = 12.5e21 */ fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16, "wlc(%g, %g, %g, %g) = %g != %g", x, T, p, L, wlc(x, T, p, L), 12.5e21); } END_TEST -@ +@ <>= TCase *tc_wlc = tcase_create("WLC"); -@ +@ <>= tcase_add_test(tc_wlc, test_wlc_at_zero); tcase_add_test(tc_wlc, test_wlc_at_half); suite_add_tcase(s, tc_wlc); -@ +@ \subsection{Model tests} @@ -1883,20 +1918,20 @@ the implementation in [[tension_balance.c]], and the unit testing in <> <> #endif /* TENSION_BALANCE_H */ -@ +@ <>= <> <> <> <> -@ +@ <>= NW_SAWSIM_MODS += tension_balance CHECK_BINS += check_tension_balance -check_tension_balance_MODS = -@ +check_tension_balance_MODS = +@ The entire force balancing problem reduces to a solving two nested one-dimensional root-finding problems. First we define the one @@ -1912,7 +1947,7 @@ double tension_balance(int num_tension_groups, void **params, double last_x, double x); -@ +@ <>= double tension_balance(int num_tension_groups, one_dim_fn_t **tension_handlers, @@ -1988,11 +2023,11 @@ double tension_balance(int num_tension_groups, F = (*tension_handlers[0])(xo, params[0]); return F; } -@ +@ <>= <> -@ +@ <>= typedef struct x_of_xo_data_struct { @@ -2001,7 +2036,7 @@ typedef struct x_of_xo_data_struct { void **handler_data; /* array of void* pointers */ double *xi; } x_of_xo_data_t; -@ +@ <>= double x_of_xo(double xo, void *pdata) @@ -2039,7 +2074,7 @@ double x_of_xo(double xo, void *pdata) #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 @@ -2052,17 +2087,17 @@ 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_relx, double min_rely, 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_relx, double min_rely, int max_steps, + double min_relx, double min_rely, int max_steps, double *pYx) { double yx, ylb, yub, x; @@ -2070,7 +2105,7 @@ double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, assert(ub >= lb); ylb = (*f)(lb, params); yub = (*f)(ub, params); - + /* check some simple cases */ if (ylb == yub) { if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */ @@ -2100,13 +2135,13 @@ double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, #endif return x; } -@ +@ The one dimensional solver needs a bracketed solution, which sometimes we don't have. Generate bracketing $x$ values through bisection or exponential growth. <>= void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub); -@ +@ <>= void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub) @@ -2141,7 +2176,7 @@ void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double } //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb); } -@ +@ \subsection{Balancing implementation} @@ -2155,7 +2190,7 @@ double length_scale = 1e-10; /* HACK */ <> <> -@ +@ <>= #include /* assert() */ @@ -2163,7 +2198,7 @@ double length_scale = 1e-10; /* HACK */ #include /* HUGE_VAL, macro constant, so don't need to link to math */ #include /* fprintf(), stdout */ #include "global.h" -@ +@ \subsection{Balancing unit tests} @@ -2172,7 +2207,7 @@ Here we check to make sure the various functions work as expected, using \citeta <> <> <
> -@ +@ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE */ @@ -2180,12 +2215,12 @@ Here we check to make sure the various functions work as expected, using \citeta <> #include "global.h" #include "tension_balance.h" -@ +@ <>= <> <> -@ +@ <>= Suite *test_suite (void) @@ -2196,7 +2231,7 @@ Suite *test_suite (void) <> return s; } -@ +@ <>= <> @@ -2216,7 +2251,7 @@ START_TEST(test_single_function) fail_unless(F = k*x, NULL); } END_TEST -@ +@ We can also test balancing two springs with different spring constants. The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}. @@ -2236,18 +2271,18 @@ START_TEST(test_double_hooke) CHECK_ERR(1e-6, Fe, F); } END_TEST -@ +@ TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above. <>= TCase *tc_tbfunc = tcase_create("tension balance function"); -@ +@ <>= tcase_add_test(tc_tbfunc, test_single_function); tcase_add_test(tc_tbfunc, test_double_hooke); suite_add_tcase(s, tc_tbfunc); -@ +@ \section{Lists} @@ -2262,7 +2297,7 @@ The interface is defined in [[list.h]], the implementation in [[list.c]], and th <> <> #endif /* LIST_H */ -@ +@ <>= <> @@ -2275,7 +2310,7 @@ The interface is defined in [[list.h]], the implementation in [[list.c]], and th <> <> <> -@ +@ <>= <> @@ -2289,13 +2324,13 @@ The interface is defined in [[list.h]], the implementation in [[list.c]], and th <> <> <> -@ +@ <>= NW_SAWSIM_MODS += list CHECK_BINS += check_list -check_list_MODS = -@ +check_list_MODS = +@ <>= typedef struct list_struct { @@ -2305,13 +2340,13 @@ typedef struct list_struct { } list_t; <> -@ +@ [[head]] and [[tail]] return pointers to the head and tail nodes of the list: <>= list_t *head(list_t *list); list_t *tail(list_t *list); -@ +@ <>= list_t *head(list_t *list) { @@ -2332,11 +2367,11 @@ list_t *tail(list_t *list) } return list; } -@ +@ <>= int list_length(list_t *list); -@ +@ <>= int list_length(list_t *list) { @@ -2351,7 +2386,7 @@ int list_length(list_t *list) } return i; } -@ +@ [[push]] inserts a node relative to the current position in the list without changing the current position. @@ -2361,10 +2396,10 @@ If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c), otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c). <>= void push(list_t **pList, void *data, int below); -@ +@ <>= void push(list_t **pList, void *data, int below) -{ +{ list_t *list, *new_node; assert(pList != NULL); list = *pList; @@ -2385,14 +2420,14 @@ void push(list_t **pList, void *data, int below) new_node->prev = list; } } -@ +@ [[pop]] removes the current domain node, moving the current position to the node after it, unless that node is [[NULL]], in which case move the current position to the node before it. <>= void *pop(list_t **pList); -@ +@ <>= void *pop(list_t **pList) { @@ -2414,14 +2449,14 @@ void *pop(list_t **pList) destroy_node(popped); return data; } -@ +@ [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data. If the cleanup function is [[NULL]], no cleanup function is called. The nodes are not popped in any particular order. <>= void list_destroy(list_t **pList, destroy_data_func_t *destroy); -@ +@ <>= void list_destroy(list_t **pList, destroy_data_func_t *destroy) { @@ -2437,12 +2472,12 @@ void list_destroy(list_t **pList, destroy_data_func_t *destroy) destroy(data); } } -@ +@ [[list_to_array]] converts a list to an array of pointers, preserving order. <>= void list_to_array(list_t *list, int *length, void ***array); -@ +@ <>= void list_to_array(list_t *list, int *pLength, void ***pArray) { @@ -2464,12 +2499,12 @@ void list_to_array(list_t *list, int *pLength, void ***pArray) *pLength = length; *pArray = array; } -@ +@ [[move]] moves the current element along the list in either direction. <>= void move(list_t **pList, int downstream); -@ +@ <>= void move(list_t **pList, int downstream) { @@ -2490,20 +2525,20 @@ void move(list_t **pList, int downstream) list = flipper; if (downstream == 0) { push(&list, data, 0); /* b>A>c>d */ - list = list->prev; /* B>a>c>d */ + list = list->prev; /* B>a>c>d */ } else { push(&list, data, 1); /* a>C>b>d */ list = list->next; /* a>c>B>d */ } *pList = list; } -@ +@ [[sort]] sorts a list from smallest to largest according to a user specified comparison. <>= typedef int (comparison_fn_t)(void *A, void *B); -@ +@ For example <>= int int_comparison(void *A, void *B) @@ -2515,10 +2550,10 @@ int int_comparison(void *A, void *B) else if (a == b) return 0; else return -1; } -@ +@ <>= void sort(list_t **pList, comparison_fn_t *comp); -@ +@ Here we use the bubble sort algorithm. <>= void sort(list_t **pList, comparison_fn_t *comp) @@ -2541,7 +2576,7 @@ void sort(list_t **pList, comparison_fn_t *comp) } *pList = list; } -@ +@ [[unique]] removes duplicates from a sorted list according to a user specified comparison. Don't do this unless you have other pointers @@ -2549,7 +2584,7 @@ any dynamically allocated data in your list, because [[unique]] just drops any non-unique data without freeing it. <>= void unique(list_t **pList, comparison_fn_t *comp); -@ +@ <>= void unique(list_t **pList, comparison_fn_t *comp) { @@ -2566,12 +2601,12 @@ void unique(list_t **pList, comparison_fn_t *comp) } *pList = list; } -@ +@ [[list_print]] prints a list to a [[FILE *]] stream. <>= void list_print(FILE *file, list_t *list, char *list_name); -@ +@ <>= void list_print(FILE *file, list_t *list, char *list_name) { @@ -2583,7 +2618,7 @@ void list_print(FILE *file, list_t *list, char *list_name) } fprintf(file, "\n"); } -@ +@ [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]]. <>= @@ -2612,14 +2647,14 @@ The user must free the data pointed to by the node on their own. <> #include "list.h" <> -@ +@ <>= #include /* assert() */ #include /* malloc(), free(), rand() */ #include /* fprintf(), stdout, FILE */ #include "global.h" /* destroy_data_func_t */ -@ +@ \subsection{List unit tests} @@ -2628,7 +2663,7 @@ Here we check to make sure the various functions work as expected, using \citeta <> <> <
> -@ +@ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE */ @@ -2636,13 +2671,13 @@ Here we check to make sure the various functions work as expected, using \citeta <> #include "global.h" #include "list.h" -@ +@ <>= <> <> <> -@ +@ <>= Suite *test_suite (void) @@ -2655,7 +2690,7 @@ Suite *test_suite (void) <> return s; } -@ +@ <>= START_TEST(test_push) @@ -2674,23 +2709,23 @@ START_TEST(test_push) } } END_TEST -@ +@ <>= TCase *tc_push = tcase_create("push"); -@ +@ <>= tcase_add_test(tc_push, test_push); suite_add_tcase(s, tc_push); -@ +@ <>= -@ +@ <>= -@ +@ <>= -@ +@ \section{Function string evaluation} @@ -2719,13 +2754,13 @@ The interface is defined in [[string_eval.h]], the implementation in [[string_ev <> <> #endif /* STRING_EVAL_H */ -@ +@ <>= NW_SAWSIM_MODS += string_eval CHECK_BINS += check_string_eval -check_string_eval_MODS = -@ +check_string_eval_MODS = +@ For an introduction to POSIX process control, see\\ \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\ @@ -2748,11 +2783,11 @@ to make it easy to change the evaluating subprocess to, say, Ruby, or the users //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL} static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}; //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL}; -@ +@ The [[i]] option lets Python know that it should run in interactive mode. In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass. In interactive mode, python acts on every instruction as soon as it is recieved. -The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply turns off the default prompting ([[ps1]] and [[ps2]]), since that is mostly for human users, and our program doesn't need it. +The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply turns off the default prompting ([[ps1]] and [[ps2]]), since that is mostly for human users, and our program doesn't need it. %The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply [[pass]]es so that the silly Python header information is not printed. We leave the prompts in, because we scan for them to determine when the output has completed. Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]]. @@ -2772,14 +2807,14 @@ We store the pipe file descriptors in an array of 4 [[int]]s, and use the follow #define STDIN 0 /* initial index of stdin pair */ #define STDOUT 2 /* initial index of stdout pair */ -@ +@ So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations. As a finishing touch, we can promote the POSIX file descriptors ([[read]]/[[write]] interface) into the more familiar [[stdio.h]] \emph{streams} ([[fprintf]]/[[fgetc]] interface) using [[fdopen]], which creates a stream from an open file descriptor. <>= extern void string_eval_setup(FILE **pIn, FILE **pOut); -@ +@ <>= void string_eval_setup(FILE **pIn, FILE **pOut) { @@ -2788,7 +2823,7 @@ void string_eval_setup(FILE **pIn, FILE **pOut) int rc; assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */ assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */ - + pid = fork(); /* split process into two copies */ if (pid == -1) { /* fork error */ @@ -2850,7 +2885,7 @@ The parent waits to confirm the child closing, recieves the child's exit status, As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits. <>= extern void string_eval_teardown(FILE **pIn, FILE **pOut); -@ +@ <>= void string_eval_teardown(FILE **pIn, FILE **pOut) @@ -2884,7 +2919,7 @@ void string_eval_teardown(FILE **pIn, FILE **pOut) } */ } -@ +@ The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals. \subsection{String evaluation implementation} @@ -2895,7 +2930,7 @@ The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals. #include "string_eval.h" <> <> -@ +@ <>= #include /* assert() */ @@ -2905,18 +2940,18 @@ The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals. #include /* pid_t */ #include /* pipe(), fork(), execvp(), alarm() */ #include /* wait() */ -@ +@ <>= <> <> -@ +@ <>= <> <> <> -@ +@ \subsection{String evaluation unit tests} @@ -2926,7 +2961,7 @@ Here we check to make sure the various functions work as expected, using \citeta <> <> <
> -@ +@ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE, atof() */ @@ -2936,12 +2971,12 @@ Here we check to make sure the various functions work as expected, using \citeta #include /* SIGKILL */ <> #include "string_eval.h" -@ +@ <>= <> <> -@ +@ <>= Suite *test_suite (void) @@ -2952,7 +2987,7 @@ Suite *test_suite (void) <> return s; } -@ +@ <>= START_TEST(test_setup_teardown) @@ -2968,7 +3003,7 @@ START_TEST(test_invalid_command) char input[80], output[80]={}; string_eval_setup(&in, &out); sprintf(input, "print ABCDefg\n"); - string_eval(in, out, input, 80, output); + string_eval(in, out, input, 80, output); string_eval_teardown(&in, &out); } END_TEST @@ -2978,7 +3013,7 @@ START_TEST(test_simple_eval) char input[80], output[80]={}; string_eval_setup(&in, &out); sprintf(input, "print 3+4*5\n"); - string_eval(in, out, input, 80, output); + string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"23\n"), NULL); string_eval_teardown(&in, &out); } @@ -2989,10 +3024,10 @@ START_TEST(test_multiple_evals) char input[80], output[80]={}; string_eval_setup(&in, &out); sprintf(input, "print 3+4*5\n"); - string_eval(in, out, input, 80, output); + string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"23\n"), NULL); sprintf(input, "print (3**2 + 4**2)**0.5\n"); - string_eval(in, out, input, 80, output); + string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"5.0\n"), NULL); string_eval_teardown(&in, &out); } @@ -3005,7 +3040,7 @@ START_TEST(test_eval_with_state) sprintf(input, "print 3+4*5\n"); fprintf(in, "A = 3\n"); sprintf(input, "print A*3\n"); - string_eval(in, out, input, 80, output); + string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"9\n"), NULL); string_eval_teardown(&in, &out); } @@ -3027,8 +3062,8 @@ suite_add_tcase(s, tc_string_eval); \section{Accelerating function evaluation} My first version-0.3 code was running very slowly. -With the options suggested in the help -([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]), +With the options suggested in the help +([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]), the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]], making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions. That's only 15 calls per solution, so the search algorithm seems reasonable. @@ -3040,13 +3075,13 @@ The number of evaluation calls could be drastically reduced, however, by impleme double accel_k(k_func_t *k, double F, environment_t *env, void *params); void free_accels(); #endif /* ACCEL_K_H */ -@ +@ <>= NW_SAWSIM_MODS += accel_k #CHECK_BINS += check_accel_k -check_accel_k_MODS = -@ +check_accel_k_MODS = +@ <>= #include /* assert() */ @@ -3090,7 +3125,7 @@ static int add_accel_k(k_func_t *k, environment_t *env, void *params) { int i=num_accels; accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels); - assert(accels != NULL); + assert(accels != NULL); accels[i].itable = interp_table_allocate(&k_wrap, &k_tol); accels[i].k = k; accels[i].env = env; @@ -3117,12 +3152,12 @@ double accel_k(k_func_t *k, double F, environment_t *env, void *params) * WARNING: we're assuming param and env are the same. */ return interp_table_eval(accels[i].itable, accels+i, F); } - } + } /* set up a new acceleration environment */ i = add_accel_k(k, env, params); return interp_table_eval(accels[i].itable, accels+i, F); } -@ +@ \section{Tension models} \label{sec.tension_models} @@ -3145,13 +3180,13 @@ The interface is defined in [[tension_model.h]], the implementation in [[tension <> <> #endif /* TENSION_MODEL_H */ -@ +@ <>= NW_SAWSIM_MODS += tension_model #CHECK_BINS += check_tension_model -check_tension_model_MODS = -@ +check_tension_model_MODS = +@ <>= TENSION_MODEL_MODS = tension_model parse list tension_balance TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \ @@ -3179,7 +3214,7 @@ For unstretchable domains. <>= {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL} -@ +@ \subsection{Constant} \label{sec.const_tension} @@ -3187,20 +3222,20 @@ For unstretchable domains. <>= <> <> -@ +@ <>= <> <> <> -@ +@ An infinitely stretchable domain providing a constant tension. <>= extern double const_tension_handler(double x, void *pdata); -@ +@ <>= double const_tension_handler(double x, void *pdata) { @@ -3222,19 +3257,19 @@ double const_tension_handler(double x, void *pdata) return F; } -@ +@ <>= typedef struct const_tension_param_struct { double F; /* tension (force) in N */ } const_tension_param_t; -@ +@ <>= extern void *string_create_const_tension_param_t(char **param_strings); extern void destroy_const_tension_param_t(void *p); -@ +@ <>= const_tension_param_t *create_const_tension_param_t(double F) @@ -3263,17 +3298,17 @@ void destroy_const_tension_param_t(void *p) extern int num_const_tension_params; extern char *const_tension_param_descriptions[]; extern char const_tension_param_string[]; -@ +@ <>= int num_const_tension_params = 1; char *const_tension_param_descriptions[] = {"tension F, N"}; 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} -@ +@ \subsection{Hooke} \label{sec.hooke} @@ -3281,14 +3316,14 @@ char const_tension_param_string[] = "0"; <>= <> <> -@ +@ <>= <> <> <> -@ +@ The tension of a single spring is given by $F=kx$ for some spring constant $k$. The behavior of a series of springs $k_i$ in series is given by @@ -3299,7 +3334,7 @@ For a simple proof, see Appendix \ref{app.math_hooke}. <>= extern double hooke_handler(double x, void *pdata); -@ +@ <>= double hooke_handler(double x, void *pdata) @@ -3327,18 +3362,18 @@ double hooke_handler(double x, void *pdata) #endif return k*x; } -@ +@ <>= typedef struct hooke_param_struct { double k; /* spring constant in N/m */ } hooke_param_t; -@ +@ <>= extern void *string_create_hooke_param_t(char **param_strings); extern void destroy_hooke_param_t(void *p); -@ +@ <>= hooke_param_t *create_hooke_param_t(double k) @@ -3365,17 +3400,17 @@ void destroy_hooke_param_t(void *p) extern int num_hooke_params; extern char *hooke_param_descriptions[]; extern char hooke_param_string[]; -@ +@ <>= int num_hooke_params = 1; char *hooke_param_descriptions[] = {"spring constant k, N/m"}; 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} -@ +@ \subsection{Worm-like chain} \label{sec.wlc} @@ -3385,14 +3420,14 @@ 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 @@ -3407,7 +3442,7 @@ $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) -{/* N m K m m */ +{/* N m K m m */ double xL=0.0; /* uses global k_B */ assert(x >= 0); assert(T > 0); assert(p > 0); assert(L > 0); @@ -3417,12 +3452,12 @@ 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); -@ +@ <>= double wlc_handler(double x, void *pdata) @@ -3450,19 +3485,19 @@ double wlc_handler(double x, void *pdata) #endif return wlc(x, data->env->T, p, L); } -@ +@ <>= typedef struct wlc_param_struct { double p; /* WLC persistence length (m) */ double L; /* WLC contour length (m) */ } wlc_param_t; -@ +@ <>= extern void *string_create_wlc_param_t(char **param_strings); extern void destroy_wlc_param_t(void *p /* wlc_param_t * */); -@ +@ <>= wlc_param_t *create_wlc_param_t(double p, double L) @@ -3502,16 +3537,16 @@ string[1] = 'x'; extern int num_wlc_params; extern char *wlc_param_descriptions[]; extern char wlc_param_string[]; -@ +@ <>= int num_wlc_params = 2; char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"}; 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} -@ +@ Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696. \subsection{Freely-jointed chain} @@ -3521,14 +3556,14 @@ Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696. <> <> <> -@ +@ <>= <> <> <> -@ +@ The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints. The tension of a chain of $N$ such links is given by @@ -3559,11 +3594,11 @@ double fjc(double x, double T, double l, int N) //printf("x %g\tL%g\tx/L %g\n", x, L, x/L); return k_B*T/l * inv_langevin(x/L); } -@ +@ <>= extern double fjc_handler(double x, void *pdata); -@ +@ <>= double fjc_handler(double x, void *pdata) { @@ -3591,19 +3626,19 @@ double fjc_handler(double x, void *pdata) #endif return fjc(x, data->env->T, l, N); } -@ +@ <>= typedef struct fjc_param_struct { double l; /* FJC link length (m) */ int N; /* FJC number of links */ } fjc_param_t; -@ +@ <>= extern void *string_create_fjc_param_t(char **param_strings); extern void destroy_fjc_param_t(void *p); -@ +@ <>= fjc_param_t *create_fjc_param_t(double l, double N) @@ -3625,23 +3660,23 @@ void destroy_fjc_param_t(void *p) if (p) free(p); } -@ +@ <>= extern int num_fjc_params; extern char *fjc_param_descriptions[]; extern char fjc_param_string[]; -@ +@ <>= int num_fjc_params = 2; char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"}; 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} -@ +@ \subsection{Tension model implementation} @@ -3653,7 +3688,7 @@ char fjc_param_string[]="0.5e-9,200"; <> <> <> -@ +@ <>= #include /* assert() */ @@ -3663,14 +3698,14 @@ char fjc_param_string[]="0.5e-9,200"; #include "global.h" #include "list.h" #include "tension_balance.h" /* oneD_*() */ -@ +@ <>= <> <> <> <> -@ +@ <>= <> @@ -3678,14 +3713,14 @@ char fjc_param_string[]="0.5e-9,200"; <> <> <> -@ +@ <>= <> <> <> <> -@ +@ \subsection{Utilities} @@ -3705,7 +3740,7 @@ range of $x$. <> <> <> -@ +@ <>= int main(int argc, char **argv) @@ -3758,7 +3793,7 @@ int main(int argc, char **argv) (*model->destructor)(params); return 0; } -@ +@ <>= #include @@ -3769,7 +3804,7 @@ int main(int argc, char **argv) #include "parse.h" #include "list.h" #include "tension_model.h" -@ +@ <>= <> @@ -3778,13 +3813,13 @@ int main(int argc, char **argv) <> <> -@ +@ <>= <> <> <> -@ +@ <>= void help(char *prog_name, @@ -3822,7 +3857,7 @@ void help(char *prog_name, printf("\t\tdefaults: %s\n", tension_models[i].params); } } -@ +@ <>= void get_options(int argc, char **argv, environment_t *env, @@ -3871,7 +3906,7 @@ void get_options(int argc, char **argv, environment_t *env, } return; } -@ +@ \section{Unfolding rate models} @@ -3881,13 +3916,13 @@ We have two main choices for determining $k$: Bell's model, which assumes the folded domains exist in a single `folded' state and make instantaneous, stochastic transitions over a free energy barrier; and Kramers' model, which treats the unfolding as a thermalized diffusion process. -We deal with these options like we dealt with the different tension models: by bundling all model-specific +We deal with these options like we dealt with the different tension models: by bundling all model-specific parameters into structures, with model specific functions for determining $k$. <>= typedef double k_func_t(double F, environment_t *env, void *params); -@ +@ Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files. The interface is defined in [[k_model.h]], the implementation in [[k_model.c]], the unit testing in [[check_k_model.c]], and some other tidbits in [[k_model_utils.c]]. @@ -3906,16 +3941,17 @@ which I don't like as much. <> <> <> +<> <> <> #endif /* K_MODEL_H */ -@ +@ <>= NW_SAWSIM_MODS += k_model CHECK_BINS += check_k_model check_k_model_MODS = parse string_eval -@ +@ [[check_k_model]] also depends on the parse module. <>= @@ -3928,7 +3964,7 @@ $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR) gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%) clean_k_model_utils : rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils -@ +@ \subsection{Null} \label{sec.null_k} @@ -3936,15 +3972,15 @@ clean_k_model_utils : For domains that are never folded. <>= -@ +@ <>= -@ +@ <>= -@ +@ <>= {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL} -@ +@ \subsection{Constant rate model} \label{sec.k_const} @@ -3958,24 +3994,24 @@ where $k_0$ is the constant unfolding rate. <>= <> <> -@ +@ <>= <> <> <> -@ +@ <>= typedef struct const_k_param_struct { double knot; /* transition rate k_0 (frac population per s) */ } const_k_param_t; -@ +@ <>= double const_k(double F, environment_t *env, void *const_k_params); -@ +@ <>= double const_k(double F, environment_t *env, void *const_k_params) { /* Returns the rate constant k in frac pop per s. */ @@ -3984,12 +4020,12 @@ double const_k(double F, environment_t *env, void *const_k_params) assert(p->knot > 0); return p->knot; } -@ +@ <>= void *string_create_const_k_param_t(char **param_strings); void destroy_const_k_param_t(void *p); -@ +@ <>= const_k_param_t *create_const_k_param_t(double knot) @@ -4026,12 +4062,12 @@ char const_k_param_string[]="1"; <>= {"const", "constant rate transitions", &const_k, num_const_k_params, const_k_param_descriptions, (char *)const_k_param_string, &string_create_const_k_param_t, &destroy_const_k_param_t} -@ +@ \subsection{Bell's model} \label{sec.bell} -Quantitatively, Bell's model gives $k$ as +Quantitatively, Bell's model gives $k$ as $$ k = k_0 \cdot e^{F dx / k_B T} \;, $$ @@ -4043,25 +4079,25 @@ $k_B T$ is the thermal energy\citep{schlierf06}. <>= <> <> -@ +@ <>= <> <> <> -@ +@ <>= typedef struct bell_param_struct { double knot; /* transition rate k_0 (frac population per s) */ double dx; /* `distance to transition state' (nm) */ } bell_param_t; -@ +@ <>= double bell_k(double F, environment_t *env, void *bell_params); -@ +@ <>= double bell_k(double F, environment_t *env, void *bell_params) { /* Returns the rate constant k in frac pop per s. @@ -4079,12 +4115,12 @@ double bell_k(double F, environment_t *env, void *bell_params) */ return p->knot * exp(F*p->dx / (k_B*env->T) ); } -@ +@ <>= void *string_create_bell_param_t(char **param_strings); void destroy_bell_param_t(void *p); -@ +@ <>= bell_param_t *create_bell_param_t(double knot, double dx) @@ -4122,10 +4158,111 @@ char bell_param_string[]="3.3e-4,0.25e-9"; <>= {"bell", "thermalized, two-state transitions", &bell_k, num_bell_params, bell_param_descriptions, (char *)bell_param_string, &string_create_bell_param_t, &destroy_bell_param_t} -@ +@ Initialized with Titin I27 parameters\citep{carrion-vazquez99b}. % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696 + +\subsection{Linker stiffness corrected Bell model} +\label{sec.kbell} + +Walton et. al showed that the Bell model constant force approximation +doesn't hold for biotin-streptavadin binding, and I extended their +results to I27 for the 2009 Biophysical Society Annual +Meeting\cite{walton08,king09}. More details to follow when I get done +with the conference... + +We adjust Bell's model to +$$ + k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;, +$$ +where $k_0$ is the unforced unfolding rate, +$F$ is the applied unfolding force, +$\kappa$ is the effective spring constant, +$dx$ is the distance to the transition state, and +$k_B T$ is the thermal energy\citep{schlierf06}. + +<>= +<> +<> +@ + +<>= +<> +<> +<> + +@ + +<>= +typedef struct kbell_param_struct { + double knot; /* transition rate k_0 (frac population per s) */ + double dx; /* `distance to transition state' (nm) */ +} kbell_param_t; +@ + +<>= +double kbell_k(double F, environment_t *env, void *kbell_params); +@ +<>= +double kbell_k(double F, environment_t *env, void *kbell_params) +{ /* Returns the rate constant k in frac pop per s. + * Takes F in N, T in K, knot in frac pop per s, and dx in m. + * uses global k_B in J/K */ + kbell_param_t *p = (kbell_param_t *) kbell_params; + assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0); + assert(p != NULL); + assert(p->knot > 0); assert(p->dx >= 0); + return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) ); +} +@ + +<>= +void *string_create_kbell_param_t(char **param_strings); +void destroy_kbell_param_t(void *p); +@ + +<>= +kbell_param_t *create_kbell_param_t(double knot, double dx) +{ + kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t)); + assert(ret != NULL); + ret->knot = knot; + ret->dx = dx; + return ret; +} + +void *string_create_kbell_param_t(char **param_strings) +{ + return create_kbell_param_t(atof(param_strings[0]), atof(param_strings[1])); +} + +void destroy_kbell_param_t(void *p /* kbell_param_t * */) +{ + if (p) + free(p); +} +@ + +<>= +extern int num_kbell_params; +extern char *kbell_param_descriptions[]; +extern char kbell_param_string[]; +@ + +<>= +int num_kbell_params = 2; +char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"}; +char kbell_param_string[]="3.3e-4,0.25e-9"; +@ + +<>= +{"kbell", "thermalized, two-state transitions, corrected for effective linker stiffness", &kbell_k, num_kbell_params, kbell_param_descriptions, (char *)kbell_param_string, &string_create_kbell_param_t, &destroy_kbell_param_t} +@ +Initialized with Titin I27 parameters\citep{carrion-vazquez99b}. +% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696 + + \subsection{Kramer's model} \label{sec.kramers} @@ -4159,14 +4296,14 @@ $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy <>= <> <> -@ +@ <>= <> <> <> -@ +@ <>= //typedef double kramers_x_func_t(double F, double T, double *E_params, int n); @@ -4187,12 +4324,12 @@ typedef struct kramers_param_struct { //double *E_params; /* parametrize the energy landscape */ //int n_E_params; /* length of E_params array */ } kramers_param_t; -@ +@ <>= extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */ extern double kramers_k(double F, environment_t *env, void *kramers_params); -@ +@ <>= double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T) { @@ -4248,12 +4385,12 @@ double kramers_k(double F, environment_t *env, void *kramers_params) //fprintf(stderr,"D = %g, lc = %g, lts = %g, Eb = %g, kbT = %g, k = %g\n", p->D, lc, lts, Eb, kbT, p->D/(lc*lts) * exp(-Eb/kbT)); return p->D/(lc*lts) * exp(-Eb/kbT); } -@ +@ <>= void *string_create_kramers_param_t(char **param_strings); void destroy_kramers_param_t(void *p); -@ +@ <>= kramers_param_t *create_kramers_param_t(double D, @@ -4313,7 +4450,7 @@ void destroy_kramers_param_t(void *p /* kramers_param_t * */) For our defaults, we'll follow \citet{schlierf06} and study ddFLN4. Schlierf and Rief used a cubic-spline interpolation routine and the double integral form of Kramers' theory, so we get to pick an actual function to approximate the energy landscape. We follow \citet{shillcock98} and use -$$ +$$ E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}}) -F x @@ -4353,7 +4490,7 @@ 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} -@ +@ \subsection{Kramer's model (natural cubic splines)} \label{sec.kramers_integ} @@ -4383,14 +4520,14 @@ We'll let the \citetalias{galassi05} natural cubic spline implementation do the <> <> <> -@ +@ <>= <> <> <> -@ +@ <>= typedef double func_t(double x, void *params); @@ -4405,7 +4542,7 @@ typedef struct kramers_integ_param_struct { double F; /* for passing into GSL evaluations */ environment_t *env; } kramers_integ_param_t; -@ +@ <>= typedef struct spline_param_struct { @@ -4415,7 +4552,7 @@ typedef struct spline_param_struct { gsl_spline *spline; /* GSL spline parameters */ gsl_interp_accel *acc; /* GSL spline acceleration data */ } spline_param_t; -@ +@ We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting. $$ @@ -4446,7 +4583,7 @@ double kramers_integ_k_integralA(double x, void *params) //fprintf(stderr, "integralA = %g\n", result); return result; } -@ +@ $$ \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv @@ -4483,7 +4620,7 @@ double kramers_integ_k_integralB(double F, environment_t *env, void *params) //fprintf(stderr, "integralB = %g\n", result); return result; } -@ +@ With the integrals taken care of, evaluating $k$ is simply $$ @@ -4491,19 +4628,19 @@ $$ $$ <>= double kramers_integ_k(double F, environment_t *env, void *kramers_params); -@ +@ <>= double kramers_integ_k(double F, environment_t *env, void *kramers_params) { /* Returns the rate constant k in frac pop per s. Takes F in N. */ kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params; return p->D/kramers_integ_k_integralB(F, env, kramers_params); } -@ +@ <>= void *string_create_kramers_integ_param_t(char **param_strings); void destroy_kramers_integ_param_t(void *p); -@ +@ <>= kramers_integ_param_t *create_kramers_integ_param_t(double D, @@ -4549,7 +4686,7 @@ Finally we have the GSL spline wrappers: <>= <> <> -@ +@ <>= double spline_eval(double x, void *spline_params) @@ -4557,7 +4694,7 @@ double spline_eval(double x, void *spline_params) spline_param_t *p = (spline_param_t *)(spline_params); return gsl_spline_eval(p->spline, x, p->acc); } -@ +@ <>= spline_param_t *create_spline_param_t(int n_knots, double *x, double *y) @@ -4612,7 +4749,7 @@ void destroy_spline_param_t(void *params) free(p); } } -@ +@ <>= extern int num_kramers_integ_params; @@ -4637,7 +4774,7 @@ char kramers_integ_param_string[]="0.2e-12,{(2e-9,0),(2.1e-9,-3.7e-20),(2.2e-9,- <>= {"kramers_integ", "thermalized, diffusion-limited transitions", &kramers_integ_k, num_kramers_integ_params, kramers_integ_param_descriptions, (char *)kramers_integ_param_string, &string_create_kramers_integ_param_t, &destroy_kramers_integ_param_t} -@ +@ Initialized with Titin I27 parameters\citep{carrion-vazquez99b}. % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696 @@ -4650,7 +4787,7 @@ Initialized with Titin I27 parameters\citep{carrion-vazquez99b}. <> <> <> -@ +@ <>= #include /* assert() */ @@ -4662,31 +4799,34 @@ Initialized with Titin I27 parameters\citep{carrion-vazquez99b}. #include #include "global.h" #include "parse.h" -@ +@ <>= <> <> +<> <> <> <> -@ +@ <>= <> <> <> +<> <> <> -@ +@ <>= <> <> <> +<> <> <> -@ +@ \subsection{Unfolding model unit tests} @@ -4697,7 +4837,7 @@ Here we check to make sure the various functions work as expected, using \citeta <> <> <
> -@ +@ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE, atof() */ @@ -4707,13 +4847,13 @@ Here we check to make sure the various functions work as expected, using \citeta <> #include "global.h" #include "k_model.h" -@ +@ <>= <> <> <> -@ +@ <>= Suite *test_suite (void) @@ -4726,19 +4866,19 @@ Suite *test_suite (void) <> return s; } -@ +@ \subsubsection{Constant} <>= TCase *tc_const_k = tcase_create("Constant k"); -@ +@ <>= tcase_add_test(tc_const_k, test_const_k_create_destroy); tcase_add_test(tc_const_k, test_const_k_over_range); suite_add_tcase(s, tc_const_k); -@ +@ <>= @@ -4750,7 +4890,7 @@ START_TEST(test_const_k_create_destroy) int i; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; - p = string_create_const_k_param_t(params); + p = string_create_const_k_param_t(params); destroy_const_k_param_t(p); } } @@ -4768,7 +4908,7 @@ START_TEST(test_const_k_over_range) env.T = T; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; - p = string_create_const_k_param_t(params); + p = string_create_const_k_param_t(params); for ( F=Fm; F>= TCase *tc_bell = tcase_create("Bell's k"); -@ +@ <>= tcase_add_test(tc_bell, test_bell_k_create_destroy); tcase_add_test(tc_bell, test_bell_k_at_zero); tcase_add_test(tc_bell, test_bell_k_at_one); suite_add_tcase(s, tc_bell); -@ +@ <>= START_TEST(test_bell_k_create_destroy) @@ -4803,7 +4943,7 @@ START_TEST(test_bell_k_create_destroy) int i; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; - p = string_create_bell_param_t(params); + p = string_create_bell_param_t(params); destroy_bell_param_t(p); } } @@ -4822,7 +4962,7 @@ START_TEST(test_bell_k_at_zero) env.T = T; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; - p = string_create_bell_param_t(params); + p = string_create_bell_param_t(params); fail_unless(bell_k(F, &env, p)==atof(knot[i]), "bell_k(%g, %g, &{%s,%s}) = %g != %s", F, T, knot[i], dx, bell_k(F, &env, p), knot[i]); @@ -4845,7 +4985,7 @@ START_TEST(test_bell_k_at_one) env.T = T; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; - p = string_create_bell_param_t(params); + p = string_create_bell_param_t(params); CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p)); //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0), // "bell_k(%g, %g, &{%s,%s}) = %g != %g", @@ -4854,16 +4994,16 @@ START_TEST(test_bell_k_at_one) } } END_TEST -@ +@ <>= -@ +@ <>= -@ +@ <>= -@ +@ <>= <> @@ -4880,7 +5020,7 @@ START_TEST(test_something) fail_unless(F = k*x, NULL); } END_TEST -@ +@ \subsection{Utilities} @@ -4897,7 +5037,7 @@ or special mode, where the operation depends on the specific model selected. <> <> <> -@ +@ <>= int main(int argc, char **argv) @@ -4968,7 +5108,7 @@ int main(int argc, char **argv) (*model->destructor)(params); return 0; } -@ +@ <>= double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x) @@ -4977,12 +5117,12 @@ double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double double E; if (k == kramers_integ_k) E = (*p->E)(x, p->E_params); - else + else E = kramers_E(0, env, model_params, x); return E; } -@ +@ <>= #include @@ -4993,7 +5133,7 @@ double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double #include "parse.h" #include "string_eval.h" #include "k_model.h" -@ +@ <>= <> @@ -5005,13 +5145,13 @@ enum mode_t {M_K_OF_F, M_SPECIAL}; <> <> -@ +@ <>= <> <> <> -@ +@ <>= void help(char *prog_name, @@ -5038,7 +5178,7 @@ void help(char *prog_name, printf(" ...\n"); printf("In special mode, the output depends on the model.\n"); printf("For models defining an energy function E(x), that function is printed\n"); - printf(" #Position (m)\tE (J)\n"); + printf(" #Distance (m)\tE (J)\n"); printf(" 0.0012\t0.0034\n"); printf(" ...\n"); printf("-m\tChange output to standard mode\n"); @@ -5057,7 +5197,7 @@ void help(char *prog_name, printf("\t\tdefaults: %s\n", k_models[i].params); } } -@ +@ <>= void get_options(int argc, char **argv, environment_t *env, @@ -5112,7 +5252,7 @@ void get_options(int argc, char **argv, environment_t *env, } return; } -@ +@ \section{\LaTeX\ documentation} @@ -5124,7 +5264,7 @@ $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) noweave -latex -index -delay $< > $@ $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR) cp -f $< $@ -@ +@ and compiled using <>= $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \ @@ -5141,14 +5281,14 @@ the remaining two [[pdflatex]] passes insert the citations and finalize the vari [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up. <>= -semi-clean_latex : +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 clean_latex : semi-clean_latex rm -f $(DOC_DIR)/sawsim.pdf -@ +@ \section{Makefile} @@ -5173,12 +5313,12 @@ SRC_DIR = src BUILD_DIR = build BIN_DIR = bin DOC_DIR = doc -# Note: Cannot use, for example, `./build', because make eats the `./' -# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...) +# Note: Cannot use, for example, `./build', because make eats the `./' +# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...) # Modules (X.c, X.h) defined in the noweb file -NW_SAWSIM_MODS = -CHECK_BINS = +NW_SAWSIM_MODS = +CHECK_BINS = # Modules defined outside the noweb file FREE_SAWSIM_MODS = interp tavl @@ -5236,8 +5376,8 @@ $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) : mkdir $@ # Copy over the source external to sawsim -# Note: Cannot use, for example, `./build', because make eats the `./' -# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...) +# Note: Cannot use, for example, `./build', because make eats the `./' +# and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...) .SECONDEXPANSION: $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\ | $(BUILD_DIR) @@ -5299,7 +5439,7 @@ $(SAWSIM_MODS:%=clean_%) : Makefile : $(SRC_DIR)/sawsim.nw notangle -Rmakefile $< | sed 's/ /\t/' > $@ -@ +@ This is fairly self-explanatory. We compile a dynamically linked version ([[sawsim]]) and a statically linked version ([[sawsim_static]]). The static version is about 9 times as big, but @@ -5324,7 +5464,7 @@ $$ x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\ x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\ F &= k_1 x_1 = k_\text{eff} x \\ - k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}} + k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}. \end{align} \end{proof}