From 99b13b65f6956c92a876ff3482fc4482daf8ed45 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Mon, 21 Jul 2008 17:30:57 +0000 Subject: [PATCH] Fixed group implementation, still need to test it. git-svn-id: svn://abax.physics.drexel.edu/sawsim/trunk@10 865a22a6-13cc-4084-8aa6-876098d8aa20 --- sawsim.nw | 1140 +++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 804 insertions(+), 336 deletions(-) diff --git a/sawsim.nw b/sawsim.nw index 3b78bcc..1c01ab7 100644 --- a/sawsim.nw +++ b/sawsim.nw @@ -73,11 +73,11 @@ \section{Introduction} -The unfolding of multi-globular protein strings is a stochastic process. -In the \prog\ program, we use Monte Carlo methods to simulate this unfolding -through an explicit time-stepping approach. -We develop a program in \lang\ to simulate probability of unfolding under -a constant extension velocity of a chain of globular domains. +The unfolding of multi-globular protein strings is a stochastic +process. In the \prog\ program, we use Monte Carlo methods to +simulate this unfolding through an explicit time-stepping approach. +We develop a program in \lang\ to simulate probability of unfolding +under a constant extension velocity of a chain of globular domains. \subsection{Related work} @@ -85,8 +85,8 @@ TODO References \subsection{About this document} -This document was written in \citetalias{sw:noweb} with code blocks in \lang\ -and comment blocks in \LaTeX. +This document was written in \citetalias{sw:noweb} with code blocks in +\lang\ and comment blocks in \LaTeX. \subsection{Change Log} @@ -96,17 +96,21 @@ and had a fixed timestep $dt$. Version 0.1 added adaptive timesteps, adjusting so that the unfolding probability for a given domain was constant. -Version 0.2 added Kramers' model unfolding rate calculations, -and a switch to select Bell's or Kramers' model. -This induced a major rewrite, introducing the tension groups abstraction, and a split to multiple \lang\ source files. -Also added a unit-test suites for sanity checking, and switched to SI units throughout. +Version 0.2 added Kramers' model unfolding rate calculations, and a +switch to select Bell's or Kramers' model. This induced a major +rewrite, introducing the tension group abstraction, and a split to +multiple \lang\ source files. Also added a unit-test suites for +sanity checking, and switched to SI units throughout. -Version 0.3 added integrated Kramers' model unfolding rate calculations. -Also replaced some of my hand-coded routines with GNU Scientific Library \citepalias{galassi05} calls. +Version 0.3 added integrated Kramers' model unfolding rate +calculations. Also replaced some of my hand-coded routines with GNU +Scientific Library \citepalias{galassi05} calls. -Version 0.4 added Python evaluation for the saddle-point Kramers' model unfolding rate calculations, -so that the model functions could be controlled from the command line. -Also added interpolating lookup tables to accelerate slow unfolding rate model evaluation. +Version 0.4 added Python evaluation for the saddle-point Kramers' +model unfolding rate calculations, so that the model functions could +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. <>= #define VERSION "0.4" @@ -147,13 +151,14 @@ Also added interpolating lookup tables to accelerate slow unfolding rate model e \section{Main} -The unfolding system is stored as a chain of \emph{domains}. -Domains can be folded, globular protein domains, unfolded protein linkers, -AFM cantilevers, or other stretchable system components. +The unfolding system is stored as a chain of \emph{domains}. Domains +can be folded, globular protein domains, unfolded protein linkers, AFM +cantilevers, or other stretchable system components. -Each domain contributes to the total tension, which depends on the chain extension. -The particular model for the tension as a function of extension is handled generally with function pointers. -So far the following models have been implemented: +Each domain contributes to the total tension, which depends on the +chain extension. The particular model for the tension as a function +of extension is handled generally with function pointers. So far the +following models have been implemented: \begin{packed_item} \item Null (Appendix \ref{sec.null_tension}), \item Constant (Appendix \ref{sec.const_tension}), @@ -162,8 +167,10 @@ So far the following models have been implemented: \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}). \end{packed_item} -The tension model and parameters used for a particular domain depend on whether the domain is folded or unfolded. -The transition rate from the folded state to the unfolded state is also handled generally with function pointers, with implementations for +The tension model and parameters used for a particular domain depend +on whether the domain is folded or unfolded. The transition rate from +the folded state to the unfolded state is also handled generally with +function pointers, with implementations for \begin{packed_item} \item Null (Appendix \ref{sec.null_k}), \item Bell's model (Appendix \ref{sec.bell}), @@ -171,12 +178,12 @@ The transition rate from the folded state to the unfolded state is also handled \item Kramers' saddle point model (Appendix \ref{sec.kramers}). \end{packed_item} -The unfolding of the chain domains is modeled in two stages. -First the tension of the chain is calculated. -Then the tension (assumed to be constant over the length of the chain), -is applied to each folded domain, and used to calculate the probability of -that subdomain unfolding in the next timestep $dt$. -Then the time is stepped forward, and the process is repeated. +The unfolding of the chain domains is modeled in two stages. First +the tension of the chain is calculated. Then the tension (assumed to +be constant over the length of the chain), is applied to each folded +domain, and used to calculate the probability of that subdomain +unfolding in the next timestep $dt$. Then the time is stepped +forward, and the process is repeated. <
>= int main(int argc, char **argv) @@ -186,18 +193,17 @@ int main(int argc, char **argv) list_t *domain_list=NULL; /* variables initialized in get_options() */ environment_t env={0}; double P, dt_max, v; - int num_folded=0; + int num_folded=0, unfolding; double x=0, dt, F; /* variables used in the simulation loop */ - one_dim_func_t *tension_handler[NUM_TENSION_GROUPS] = {0}; - get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_GROUPS, + get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded); - setup(tension_handler); + setup(); if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n"); if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n"); while (num_folded > 0) { - F = find_tension(tension_handler, domain_list, &env, x); - dt = determine_dt(tension_handler, domain_list, &env, x, P, dt_max, v); - random_unfoldings(domain_list, F, dt, &env, &num_folded); + F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding); + dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v); + unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded); if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F); x += v*dt; } @@ -211,9 +217,9 @@ int main(int argc, char **argv) [[determine_dt]] in Section \ref{sec.adaptive_dt}, and [[random_unfoldings]] in Section \ref{sec.unfolding_rate}. -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. +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; @@ -240,84 +246,157 @@ typedef struct environment_struct { \subsection{Tension} \label{sec.find_tension} -Because the stretched system may be made up of several parts (folded domains, unfolded domains, spring-like cantilever, \ldots), -we will assign the domains to groups. For example, a domain might be in the [[NULL]] group when it's folded and in the worm-like chain group when it is unfolded. -The domains are assumed to be commutative, so ordering is ignored. -The interactions between the groups are assumed to be linear, but the interactions between domains of the same group need not be. -This allows for non-linear group models such as th worm-like or freely-jointed chains. -Each handler function receives a list of domains matching it's group number. +Because the stretched system may be made up of several parts (folded +domains, unfolded domains, spring-like cantilever, \ldots), we will +assign the domains to models and groups. The models are used to +determine a domain's tension (Hookean spring, WLC, \ldots). Groups +will represent collections of domains which the model should treat as +a single entity. A domain may belong to different groups or models +depending on its state. For example, a domain might be modeled as a +freely-jointed chain when it iss folded and as a worm-like chain class +when it is unfolded. The domains are assumed to be commutative, so +ordering is ignored. The interactions between the groups are assumed +to be linear, but the interactions between domains of the same group +need not be. This allows for non-linear group models such as th +worm-like or freely-jointed chains. Each model has a tension handler +function, which gives the tension $F$ needed to stretch a given group +of domains a distance $x$. + +To understand the purpose of groups, consider a chain of proteins +where two folded proteins $A$ and $B$ are modeled as WLCs with +persistence length $p_f$ and two unfolded proteins $C$ and $D$ are +modeled as WLCs with persistence length $p_u$. The proteins are +attached to a cantilever $E$ qof spring constant $κ$. The string +would get split into three lists: +\begin{center} + \begin{tabular}{llll} + Model & Group & List & Tension \\ + WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B}$ \\ + WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D}$ \\ + Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ}$ \\ + \end{tabular} +\end{center} +Note that group numbers only matter \emph{within} model classes, since +grouping domains with seperate models doesn't make sense. + <>= -enum tension_group_t {TG_NULL=0, - TG_CONST, - TG_HOOKE, - TG_WLC, - TG_FJC, -}; -#define NUM_TENSION_GROUPS 5 +#define NUM_TENSION_MODELS 5 + +@ + +<>= +extern one_dim_fn_t *tension_handlers[]; +@ + +<>= +one_dim_fn_t *tension_handlers[] = { + NULL, + &const_tension_handler, + &hooke_handler, + &wlc_handler, + &fjc_handler, + }; + +@ +<>= +<> @ <>= typedef struct tension_handler_data_struct { - list_t *group; + /* int num_domains; */ + /* domain_t *domains; */ + list_t *group; environment_t *env; - void *persist; + 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 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\; \forall i,j$. -For the moment, we will restrict our group boundaries to a single dimension, so $\sum_i x_i = x$, our total extension (it is this restriction that causes the groups to interact linearly). -We'll also restrict our extensions to all be positive. -With these restrictions, the problem of balancing the tensions reduces to solving a system of $N+1$ possibly non-linear equations with $N$ unknowns, -where $N$ is the number of active groups. -In general this can be fairly complicated, but without much loss of practicality we can restrict ourselves to strictly monotonically increasing, non-negative tension functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much simpler. -For example, it guarantees the existence of a unique, real solution for finite forces. +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 +each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\; +\forall i,j$. Note that there may be several groups within each model +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]]. <>= -double find_tension(one_dim_func_t **tension_handler, list_t *domain_list, environment_t *env, double x) +double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding) { - static tension_handler_data_t data[NUM_TENSION_GROUPS] = {0}; - static void *pdata[NUM_TENSION_GROUPS] = {0}; - static double xi[NUM_TENSION_GROUPS] = {0}, last_x = 0; - static int active_groups[NUM_TENSION_GROUPS] = {0}; - int i, new_active_groups=0; + static int num_active_groups=0; + static one_dim_fn_t **per_group_handlers = NULL; + static void **per_group_params = NULL; + static double last_x; + tension_handler_data_t *tdata; double F; + int i; - for (i=0; igroup, " parameters"); +#endif + tdata->env = env; + /* tdata->persist continues unchanged */ + } + last_x = -1.0; + } /* else roll with the current statics */ + + F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x); last_x = x; return F; } -@ -See Appendix \ref{app.tension_balance} for the tension balancing implementation. - -Each group must define a parameter structure if the tension model is parametrized, +@ For the moment, we will restrict our group boundaries to a single +dimension, so $\sum_i x_i = x$, our total extension (it is this +restriction that causes the groups to interact linearly). We'll also +restrict our extensions to all be positive. With these restrictions, +the problem of balancing the tensions reduces to solving a system of +$N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is +the number of active groups. In general this can be fairly +complicated, but without much loss of practicality we can restrict +ourselves to strictly monotonically increasing, non-negative tension +functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much +simpler. For example, it guarantees the existence of a unique, real +solution for finite forces. See Appendix \ref{app.tension_balance} +for the tension balancing implementation. + +Each group must define a parameter structure if the tension model is +parametrized, a function to create the parameter structure, a function to destroy the parameter structure, and - a tension function of type [[one_dim_func_t]] that receives a [[tension_handler_data_t *]] as it's data argument. -The tension-specific parameter structures are contained in the domain group field. -See the Section \ref{app.model_selection} for a discussion on the command line framework. -See the worm-like chain implementation for an example (Section \ref{sec.wlc}). + a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument. +The tension-specific parameter structures are contained in the domain +group field. See the Section \ref{app.model_selection} for a +discussion on the command line framework. See the worm-like chain +implementation for an example (Section \ref{sec.wlc}). -The major limitation of our approach is that all of our tension models are -Markovian (memory-less), which is not really a problem for the chains -(see \citet{evans99} for WLC thermalization rate calculations). +The major limitation of our approach is that all of our tension models +are Markovian (memory-less), which is not really a problem for the +chains (see \citet{evans99} for WLC thermalization rate calculations). \subsection{Unfolding rate} \label{sec.unfolding_rate} @@ -331,8 +410,8 @@ With the rate constant $k$ defined by $$ \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}]. $$ -So the probability of a given protein unfolding in a short time $dt$ is given - by +So the probability of a given protein unfolding in a short time $dt$ +is given by $$ P = k \cdot dt. $$ @@ -348,24 +427,21 @@ int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain) } @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}. -Only one domain can unfold in each timestep, -because the timescale of a domain unfolding $dt_u$ -is assumed to be less than the simulation timestep $dt$, -so a domain will completely unfold in a single timestep. -We adapt our timesteps to keep the probability of a single domain unfolding -low, so the probability of two domains unfolding in the same timestep is -negligible. -This approach breaks down as the adaptive timestep scheme approaches -$dt \sim dt_u$, but $dt_u \sim 1\U{$\mu$s}$ for Ig27-like proteins -\citep{klimov00}, so this shouldn't be a problem. -To reassure yourself, you can ask the simulator to print the smallest timestep -that was used. - +Only one domain can unfold in each timestep, because the timescale of +a domain unfolding $dt_u$ is assumed to be less than the simulation +timestep $dt$, so a domain will completely unfold in a single +timestep. We adapt our timesteps to keep the probability of a single +domain unfolding low, so the probability of two domains unfolding in +the same timestep is negligible. This approach breaks down as the +adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim +1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this +shouldn't be a problem. To reassure yourself, you can ask the +simulator to print the smallest timestep that was used with TODO. <>= -void random_unfoldings(list_t *domain_list, double tension, - double dt, environment_t *env, - int *num_folded_domains) -{ +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 && domain_unfolds(tension, dt, env, D(domain_list))) { @@ -373,37 +449,40 @@ void random_unfoldings(list_t *domain_list, double tension, fprintf(stdout, "%g\n", tension); D(domain_list)->state = DS_UNFOLDED; (*num_folded_domains)--; - break; /* our one domain has unfolded, stop looking for others */ + return 1; /* our one domain has unfolded, stop looking for others */ } domain_list = domain_list->next; } + return 0; } @ \subsection{Adaptive timesteps} \label{sec.adaptive_dt} -We'd like to pick $dt$ so the probability of unfolding in the next timestep is -small. If we don't adapt our timestep, we risk breaking our approximation -that $F(x' \in [x, x+v \cdot dt]) = F(x)$ -or that only one domain unfolds in a given timestep. -Because $F(x)$ is monotonically increasing, excessively large timesteps will -lead to erroneously large unfolding forces. -The simulation would have been accurate for sufficiently small fixed timesteps, -but adaptive timesteps will allow us to move through low tension regions in -fewer steps, leading to a more efficient simulation. - -The actual adaptive timestep implementation is not particularly interesting, -since we are only required to reduce $dt$ to somewhere below a set threshold, -so I've removed it to Appendix \ref{app.adaptive_dt}. +We'd like to pick $dt$ so the probability of unfolding in the next +timestep is small. If we don't adapt our timestep, we risk breaking +our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that +only one domain unfolds in a given timestep. Because $F(x)$ is +monotonically increasing, excessively large timesteps will lead to +erroneously large unfolding forces. The simulation would have been +accurate for sufficiently small fixed timesteps, but adaptive +timesteps will allow us to move through low tension regions in fewer +steps, leading to a more efficient simulation. + +The actual adaptive timestep implementation is not particularly +interesting, since we are only required to reduce $dt$ to somewhere +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 allows interchangeable models, -and we are currently focusing on the simulation itself, -so we remove the actual model implementation to the appendices. + +The models provide the physics of the simulation, but the simulation +allows interchangeable models, and we are currently focusing on the +simulation itself, so we remove the actual model implementation to the +appendices. The tension models are defined in Appendix \ref{sec.tension_models} and unfolding models are defined in Appendix \ref{sec.k_models}. @@ -426,8 +505,10 @@ and unfolding models are defined in Appendix \ref{sec.k_models}. \subsection{Model selection} \label{app.model_selection} -The main difficulty with the command line interface in \prog\ is developing an intuitive method for accessing the possibly large number of available models. -We'll treat this generally by defining an array of available models, containing their important parameters. +The main difficulty with the command line interface in \prog\ is +developing an intuitive method for accessing the possibly large number +of available models. We'll treat this generally by defining an array +of available models, containing their important parameters. <>= <> @@ -446,21 +527,24 @@ typedef void destroy_data_func_t(void *param_struct); \subsubsection{Tension} +To access the various tension models from the command line, we define +a structure that contains all the neccessary information about the +model. <>= typedef struct tension_model_getopt_struct { - enum tension_group_t tg_group; - char *name; - char *description; - int num_params; - char **param_descriptions; - char *params; - create_data_func_t *creator; - destroy_data_func_t *destructor; + one_dim_fn_t *handler; /* fn implementing the model on a 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 */ + char **param_descriptions; /* descriptions of extra params */ + char *params; /* default values for the extra params */ + 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_GROUPS] = { +tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = { <>, <>, <>, @@ -520,11 +604,11 @@ void help(char *prog_name, double P, double dt_max, double v, printf("Model options:\n"); printf("The domains exist in either the folded or unfolded state\n"); printf("The following options change the default behavior in each state.\n"); - printf("-m model\tFolded tension model (currently %s)\n", + printf("-m model[,group]\tFolded tension model (currently %s)\n", tension_models[folded_tension_model].name); printf("-a args\tFolded tension model argument string (currently %s)\n", tension_models[folded_tension_model].params); - printf("-M model\tUnfolded tension model (currently %s)\n", + printf("-M model[,group]\tUnfolded tension model (currently %s)\n", tension_models[unfolded_tension_model].name); printf("-A args\tUnfolded tension model argument string (currently %s)\n", tension_models[unfolded_tension_model].params); @@ -564,8 +648,11 @@ void help(char *prog_name, double P, double dt_max, double v, printf("\t\t\t%s\n", k_models[i].param_descriptions[j]); printf("\t\tdefaults: %s\n", k_models[i].params); } - printf("Example usage: (1 spring and 8 folded domains)\n"); + printf("Examples:\n"); + printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n"); printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name); + 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); } @ @@ -582,12 +669,23 @@ void get_options(int argc, char **argv, char *prog_name = NULL; char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh"; int ftension_model=0, utension_model=0, k_model=0; - int i, n; + char *ftension_params=NULL, *utension_params=NULL; + int i, n, ftension_group=0, utension_group=0; + int num_strings; + char **strings; extern char *optarg; extern int optind, optopt, opterr; assert (argc > 0); +#ifdef DEBUG + for (i=0; iT = atof(optarg); break; case 'C': env->T = atof(optarg)+273.15; break; case 'm': - ftension_model = index_tension_model(n_tension_models, tension_models, optarg); + parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings); + assert(num_strings == 1 || num_strings == 2); + if (num_strings == 1) + ftension_group = 0; + else + ftension_group = atoi(strings[1]); + ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]); + ftension_params = tension_models[ftension_model].params; + free_parsed_list(num_strings, strings); break; case 'a': - tension_models[ftension_model].params = optarg; + ftension_params = optarg; break; case 'M': - utension_model = index_tension_model(n_tension_models, tension_models, optarg); + parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings); + assert(num_strings == 1 || num_strings == 2); + if (num_strings == 1) + utension_group = 0; + else + utension_group = atoi(strings[1]); + utension_model = index_tension_model(n_tension_models, tension_models, strings[0]); + utension_params = tension_models[utension_model].params; + free_parsed_list(num_strings, strings); break; case 'A': - tension_models[utension_model].params = optarg; + utension_params = optarg; break; case 'k': k_model = index_k_model(n_k_models, k_models, optarg); @@ -625,7 +739,14 @@ void get_options(int argc, char **argv, n = atoi(optarg); assert(n > 0); for (i=0; i 0); for (i=0; iparams, SEP, '{', '}', \ - &num_param_args, ¶m_args); \ - if (num_param_args != model->num_params) { \ - fprintf(stderr, \ - "%s model %s expected %d params," \ - role, model->name, model->num_params); \ - fprintf(stderr, \ - "not the %d params in '%s'\n", \ - num_param_args, model->params); \ - assert(num_param_args == model->num_params); \ - } \ - if (model->creator) \ - param_pointer = (*model->creator)(param_args); \ - else \ - param_pointer = NULL; \ - free_parsed_list(num_param_args, param_args); \ +#define INIT_MODEL(role, model, param_string, param_pointer) \ + do { \ + parse_list_string(param_string, SEP, '{', '}', \ + &num_param_args, ¶m_args); \ + if (num_param_args != model->num_params) { \ + fprintf(stderr, \ + "%s model %s expected %d params," \ + role, model->name, model->num_params); \ + fprintf(stderr, \ + "not the %d params in '%s'\n", \ + num_param_args, param_string); \ + assert(num_param_args == model->num_params); \ + } \ + if (model->creator) \ + param_pointer = (*model->creator)(param_args); \ + else \ + param_pointer = NULL; \ + free_parsed_list(num_param_args, param_args); \ } while (0); @ <>= void *generate_domain(enum domain_state_t state, tension_model_getopt_t *folded_model, + int folded_group, + char *folded_param_string, tension_model_getopt_t *unfolded_model, + int unfolded_group, + char *unfolded_param_string, k_model_getopt_t *k_model) { - void *ftension_params, *utension_params, *k_params; + void *folded_params, *unfolded_params, *k_params; int num_param_args; /* for INIT_MODEL() */ char **param_args; /* for INIT_MODEL() */ #ifdef DEBUG fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded"); - fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\", \"%s\", u: \"%s\", \"%s\")\n", + fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n", k_model->name, k_model->params, - folded_model->name, folded_model->params, - unfolded_model->name, unfolded_model->params); + folded_model->name, folded_model->handler, folded_group, folded_param_string, + unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string); #endif - INIT_MODEL("folded", folded_model, ftension_params); - INIT_MODEL("unfolded", unfolded_model, utension_params); - INIT_MODEL("k", k_model, k_params); + INIT_MODEL("folded", folded_model, folded_param_string, folded_params); + INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params); + INIT_MODEL("k", k_model, k_model->params, k_params); return create_domain(state, - k_model->k, k_params, - k_model->destructor, - folded_model->tg_group, ftension_params, - folded_model->destructor, - unfolded_model->tg_group, utension_params, - unfolded_model->destructor); + k_model->k, k_params, k_model->destructor, + folded_model->handler, folded_group, folded_params, folded_model->destructor, + unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor + ); } @ @@ -809,14 +939,9 @@ We include [[math.h]], so don't forget to link to the libm with `-lm'. @ <>= -void setup(one_dim_func_t **tension_handler) +void setup(void) { srand(getpid()*time(NULL)); /* seed rand() */ - tension_handler[TG_NULL] = NULL; - tension_handler[TG_CONST] = &const_tension_handler; - tension_handler[TG_HOOKE] = &hooke_handler; - tension_handler[TG_WLC] = &wlc_handler; - tension_handler[TG_FJC] = &fjc_handler; } @ @@ -886,7 +1011,7 @@ and we expect the next step to be slightly shorter than the previous one. 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(one_dim_func_t **tension_handler, +double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, double target_prob, double max_dt, double v) @@ -899,7 +1024,7 @@ double search_dt(one_dim_func_t **tension_handler, double F, k, dtCur, dtU, dtUCur, dtL, dt; /* get upper bound using the current position */ - F = find_tension(tension_handler, domain_list, env, x); /* BUG. repeated calculation */ + F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */ //printf("Start with x = %g (F = %g)\n", x, F); k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params); //printf("x %g\tF %g\tk %g\n", x, F, k); @@ -914,11 +1039,11 @@ double search_dt(one_dim_func_t **tension_handler, /* The dt determined above may produce illegitimate forces or ks. * Reduce the upper bound until we have valid ks. */ dt = dtU; - F = find_tension(tension_handler, domain_list, env, x+v*dt); + F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0); while (F == HUGE_VAL) { /* reduce step until we hit a valid force */ dtU /= 2.0; dt = dtU; - F = find_tension(tension_handler, domain_list, env, x+v*dt); + 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); @@ -927,7 +1052,7 @@ double search_dt(one_dim_func_t **tension_handler, 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(tension_handler, domain_list, env, x+v*dt); + 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); @@ -940,7 +1065,7 @@ double search_dt(one_dim_func_t **tension_handler, /* dtU wasn't safe, lets see what would be. */ while (dtU > 1.1*dtL) { /* until the range is reasonably small */ dt = (dtU + dtL) / 2.0; - F = find_tension(tension_handler, domain_list, env, x+v*dt); + F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0); //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU); k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params); dtCur = target_prob / k; @@ -961,7 +1086,7 @@ 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. <>= <> -double determine_dt(one_dim_func_t **tension_handler, list_t *domain_list, +double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, double target_prob, double dt_max, double v) { /* Returns the timestep dt in seconds. @@ -975,7 +1100,7 @@ double determine_dt(one_dim_func_t **tension_handler, list_t *domain_list, //return 0.5e-9/v; while (domain_list != NULL) { if (D(domain_list)->state == DS_FOLDED) { - new_dt = search_dt(tension_handler, domain_list, env, x, target_prob, dt, v); + new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v); dt = MIN(dt, new_dt); } domain_list = domain_list->next; @@ -986,12 +1111,13 @@ double determine_dt(one_dim_func_t **tension_handler, list_t *domain_list, \subsection{Domain data} -Currently domains exist in two states, folded and unfolded, -and the only allowed transitions are folded $\rightarrow$ unfolded. -Of course, it wouldn't be too complicated to extent this to a multi-state system, -with an array containing the domains group for each possible state, -and a matrix of transition-rate-calculating functions. -However, at this point such generality seems unnecessary at this point. +Currently domains exist in two states, folded and unfolded, and the +only allowed transitions are folded $\rightarrow$ unfolded. Of +course, it wouldn't be too complicated to extent this to a multi-state +system, with an array containing the domains group for each possible +state, and a matrix of transition-rate-calculating functions. +However, at this point such generality seems unnecessary at this +point. <>= enum domain_state_t {DS_FOLDED, DS_UNFOLDED @@ -999,9 +1125,11 @@ enum domain_state_t {DS_FOLDED, typedef struct domain_struct { enum domain_state_t state; - enum tension_group_t folded_group; - enum tension_group_t unfolded_group; - k_func_t *k; /* function returning unfolding rate */ + one_dim_fn_t *folded_handler; + int folded_group; + one_dim_fn_t *unfolded_handler; + int unfolded_group; + k_func_t *k; /* function returning unfolding rate */ void *folded_params; /* pointer to folded parameters */ void *unfolded_params; /* pointer to unfolded parameters */ void *k_params; /* pointer to k parameters */ @@ -1029,10 +1157,12 @@ domain_t *create_domain(enum domain_state_t state, k_func_t *k, void *k_params, destroy_data_func_t *destroy_k, - enum tension_group_t folded_group, + one_dim_fn_t *folded_handler, + int folded_group, void *folded_params, destroy_data_func_t *destroy_folded, - enum tension_group_t unfolded_group, + one_dim_fn_t *unfolded_handler, + int unfolded_group, void *unfolded_params, destroy_data_func_t *destroy_unfolded) { @@ -1041,17 +1171,23 @@ domain_t *create_domain(enum domain_state_t state, if (state == DS_FOLDED) { assert(k != NULL); /* the pointer points somewhere valid */ assert(*k != NULL); /* and there is something useful there */ - } + } else + assert(state == DS_UNFOLDED); ret->state = state; + ret->folded_handler = folded_handler; ret->folded_group = folded_group; + ret->unfolded_handler = unfolded_handler; ret->unfolded_group = unfolded_group; ret->k = k; - ret->k_params = k_params; - ret->destroy_k = destroy_k; ret->folded_params = folded_params; ret->unfolded_params = unfolded_params; + ret->k_params = k_params; ret->destroy_folded = destroy_folded; ret->destroy_unfolded = destroy_unfolded; + ret->destroy_k = destroy_k; +#ifdef DEBUG + fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state); +#endif return ret; } @@ -1079,41 +1215,182 @@ void destroy_domain_list(list_t *domain_list) } @ -\subsection{Group handling} +\subsection{Domain model and group handling} <>= +<> <> +<> +<> +<> <> +<> +@ + +<>= +one_dim_fn_t *get_tension_handler(domain_t *domain) +{ + if (domain->state == DS_FOLDED) + return domain->folded_handler; + else { + if (domain->state != DS_UNFOLDED) + fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state); + assert(domain->state == DS_UNFOLDED); + return domain->unfolded_handler; + } +} @ <>= -enum tension_group_t get_group(domain_t *domain) +int get_group(domain_t *domain) { if (domain->state == DS_FOLDED) return domain->folded_group; else { + if (domain->state != DS_UNFOLDED) + fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state); assert(domain->state == DS_UNFOLDED); 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 +groups. We define a function [[find_possible_groups]] to search for +possible (not neccessarily active) groups. It can search for +subgroups of a single [[handler]], or by setting [[handler]] to +[[NULL]], subgroups of \emph{any} handler. It returns a list of group +integers. +<>= +list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) { + list_t *ret = NULL; + list = head(list); + while (list != NULL) { + if (handler == NULL || get_tension_handler(D(list)) == handler) { + push(&ret, &D(list)->folded_group, 1); + push(&ret, &D(list)->unfolded_group, 1); + } + list = list->next; + } + + if (ret == NULL) + return ret; /* no groups with this handler, no processing needed */ + + /* sort the ret list, and remove duplicates */ + sort(&ret, &int_comparison); + unique(&ret, &int_comparison); +#ifdef DEBUG + fprintf(stderr, "handler %p possible groups:", handler); + list = head(ret); + while (list != NULL) { + fprintf(stderr, " %d", *((int *)list->d)); + list = list->next; + } + fprintf(stderr, "\n"); +#endif + return ret; +} +@ + +Search a [[list]] of domains, and determine whether a particular model +class and group number combination is active. +<>= +int is_group_active(list_t *list, one_dim_fn_t *handler, int group) +{ + list = head(list); + while (list != NULL) { + if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) + return 1; + list = list->next; + } + return 0; +} +@ + +Search a [[list]] of domains, and return all domains matching a +particular model class and group number. <>= -list_t *get_group_list(list_t *list, enum tension_group_t group) +list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group) { list_t *ret = NULL; list = head(list); while (list != NULL) { - if (get_group(D(list)) == group) - push(&ret, D_TP(list)); /* add a pointer to the appropriate tension parameters to our new list. */ + if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) { + push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */ +#ifdef DEBUG + fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret); +#endif + } list = list->next; } 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 and node data popped off when destroying the group lists. -It will all get cleaned up when the main domain list is destroyed. +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 +up when the main domain list is destroyed. + +Put all this together to scan through a list of domains and construct +an array of all the active groups. +<>= +void get_active_groups(list_t *list, + int num_tension_handlers, one_dim_fn_t **tension_handlers, + int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups) +{ + void **active_groups=NULL; + one_dim_fn_t *handler, **per_group_handlers=NULL; + int i, num_active_groups, *group; + list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL; + + /* get all the active groups in a list */ + for (i=0; id); + list_print(stderr, active_group_list, "active group"); +#endif + } + } + } + + /* 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); + per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups); + assert(per_group_handlers != NULL); + + /* populate the array we'll be returning */ + active_groups_list = head(active_groups_list); + for (i=0; igroup = active_group_list; + ((tension_handler_data_t *)active_groups[i])->env = NULL; + ((tension_handler_data_t *)active_groups[i])->persist = NULL; + } + assert(active_groups_list == NULL); + assert(per_group_handlers_list == NULL); + + *pNum_active_groups = num_active_groups; + *pPer_group_handlers = per_group_handlers; + *pActive_groups = active_groups; +} +@ \section{String parsing} @@ -1153,15 +1430,18 @@ extern void parse_list_string(char *string, char sep, char deeper, char shallowe 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. +[[parse_list_string]] allocates memory, don't forget to free it +afterward with [[free_parsed_list]]. It does not alter the original. -The string may start off with a [[deeper]] character (i.e. [["{x,y}"]]), -and when it does, brace stripping will set leave [[x,y]], where the pointer is one character in on the copied string. -However, when we go to free the memory, we need a pointer to the beginning of the string. -In order to accommodate this for a string with $N$ argument, allocate a pointer array with $N+1$ elements, -let the first $N$ elements point to the separated arguments, -and let the last element point to the start of the copied string regardless of braces. +The string may start off with a [[deeper]] character +(i.e. [["{x,y}"]]), and when it does, brace stripping will set leave +[[x,y]], where the pointer is one character in on the copied string. +However, when we go to free the memory, we need a pointer to the +beginning of the string. In order to accommodate this for a string +with $N$ argument, allocate a pointer array with $N+1$ elements, let +the first $N$ elements point to the separated arguments, and let the +last element point to the start of the copied string regardless of +braces. <>= /* TODO, split out into parse.hc */ static int next_delim_index(char *string, char sep, char deeper, char shallower) @@ -1530,9 +1810,14 @@ suite_add_tcase(s, tc_determine_dt); \section{Balancing group extensions} \label{app.tension_balance} -For a given total extension $x$ (of the piezo), the various domain groups (WLC, FJC, Hookean springs, \ldots) extend different amounts, and we need to tweak the portion that each extends to balance the tension amongst the active groups. -Since the tension balancing is separable from the bulk of the simulation, we break it out into a separate module. -The interface is defined in [[tension_balance.h]], the implementation in [[tension_balance.c]], and the unit testing in [[check_tension_balance.c]] +For a given total extension $x$ (of the piezo), the various domain +models (WLC, FJC, Hookean springs, \ldots) and groups extend different +amounts, and we need to tweak the portion that each extends to balance +the tension amongst the active groups. Since the tension balancing is +separable from the bulk of the simulation, we break it out into a +separate module. The interface is defined in [[tension_balance.h]], +the implementation in [[tension_balance.c]], and the unit testing in +[[check_tension_balance.c]] <>= <> @@ -1559,100 +1844,88 @@ clean_tension : rm -f tension_balance.c tension_balance.h check_tension_balance.c check_tension_balance @ -The entire force balancing problem reduces to a solving two nested one-dimensional root-finding problems. -First we define the one dimensional tension $F(x_0)$ (where $i=0$ is the index of the first non-empty group). -$x(x_0)$ is also strictly monotonically increasing, so we can solve for $x_0$ such that $\sum_i x_i = x$. +The entire force balancing problem reduces to a solving two nested +one-dimensional root-finding problems. First we define the one +dimensional tension $F(x_0)$ (where $i=0$ is the index of the first +non-empty group). $x(x_0)$ is also strictly monotonically increasing, +so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]] +stores the last successful value of $x$, since we expect to be taking +small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]] +indicates that the groups have changed. <>= double tension_balance(int num_tension_groups, - one_dim_func_t **tension_handler, - void **params, int *active, - int new_active_groups, double *xi, - double last_x, double x); + one_dim_fn_t **tension_handlers, + void **params, + double last_x, + double x); @ <>= double tension_balance(int num_tension_groups, - one_dim_func_t **tension_handler, - void **params, int *active, - int new_active_groups, double *xi, - double last_x, double x) -{ /* xi initialized to x values for groups at last x, - * returned as x values for groups at current x */ - double F, xo; - one_dim_func_t **active_handlers=NULL; - void **active_params=NULL; - double *active_xi=NULL; - int i, active_groups=0; - x_of_xo_data_t x_xo_data; + one_dim_fn_t **tension_handlers, + void **params, + double last_x, + double x) +{ + static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL}; + double F, x_guess, xo, lb, ub; double min_dx=1e-10, min_dy=1e-15; - int max_steps=100; - double lb, ub; + int max_steps=100, i; assert(num_tension_groups > 0); - active_handlers = (one_dim_func_t **)malloc(sizeof(one_dim_func_t *)*num_tension_groups); - assert(active_handlers != NULL); - active_params = (void **)malloc(sizeof(void *)*num_tension_groups); - assert(active_params != NULL); - active_xi = (double *)malloc(sizeof(double)*num_tension_groups); - assert(active_xi != NULL); - - for (i=0; i 0); + + if (num_tension_groups == 1) { /* only one group, no balancing required */ xo = x; - active_xi[0] = x; } else { //printf("balancing for x = %g with ", x); - //for(i=0; i last_x) { - lb = active_xi[0]; - ub = active_xi[0]+ x-last_x; /* apply all change to x0 */ + lb = x_xo_data.xi[0]; + ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */ } else if (x < last_x) { - lb = active_xi[0]- (x-last_x); /* apply all change to x0 */ - ub = active_xi[0]; + lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */ + ub = x_xo_data.xi[0]; } else { /* x == last_x */ printf("not moving\n"); - lb= active_xi[0]; - ub= active_xi[0]; + lb= x_xo_data.xi[0]; + ub= x_xo_data.xi[0]; } } //printf("lb %g,\tub %g\n", lb, ub); xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_dx, min_dy, max_steps, NULL); - } - F = (*active_handlers[0])(xo, active_params[0]); - /* go back through and place the active xi data in the complete xi array */ - active_groups = 0; - for (i=0; i 1 && 0) { + printf("balanced for x = %g with ", x); + for(i=0; i 1 && 0) { - printf("balanced for x = %g with ", x); - for(i=0; i>= typedef struct x_of_xo_data_struct { int n_groups; - one_dim_func_t **tension_handler; /* array of fn pointers */ - void **handler_data; /* array of void* pointers */ + one_dim_fn_t **tension_handler; /* array of fn pointers */ + void **handler_data; /* array of void* pointers */ double *xi; } x_of_xo_data_t; @ @@ -1674,7 +1947,7 @@ typedef struct x_of_xo_data_struct { double x_of_xo(double xo, void *pdata) { x_of_xo_data_t *data = (x_of_xo_data_t *)pdata; - double F, x=0, xi, lb, ub; + double F, x=0, xi, xi_guess, lb, ub; int i; double min_dx=1e-10, min_dy=1e-14; int max_steps=100; @@ -1684,7 +1957,9 @@ 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++) { - oneD_bracket(data->tension_handler[i], data->handler_data[i], F, data->xi[i], &lb, &ub); + if (data->xi[i] == -1) xi_guess = 1.0; + else xi_guess = data->xi[i]; + oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub); xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_dx, min_dy, max_steps, NULL); data->xi[i] = xi; x += xi; @@ -1698,17 +1973,17 @@ Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited number of ste Simple bisection, so it's robust and converges fairly quickly. <>= /* equivalent to gsl_func_t */ -typedef double one_dim_func_t(double x, void *params); +typedef double one_dim_fn_t(double x, void *params); @ <>= -double oneD_solve(one_dim_func_t *f, void *params, double y, double lb, double ub, +double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, double min_dx, double min_dy, int max_steps, double *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_func_t *f, void *params, double y, double lb, double ub, +double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, double min_dx, double min_dy, int max_steps, double *pYx) { @@ -1744,11 +2019,11 @@ double oneD_solve(one_dim_func_t *f, void *params, double y, double lb, double u 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_func_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); @ <>= -void oneD_bracket(one_dim_func_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) { double yx, step, x=xguess; int i=0; @@ -1835,7 +2110,7 @@ double hooke(void *pK, double x) START_TEST(test_single_function) { double k=5, x=3, last_x=2.0, F; - one_dim_func_t *handlers[] = {&hooke}; + one_dim_fn_t *handlers[] = {&hooke}; void *data[] = {&k}; double xi[] = {0.0}; int active[] = {1}; @@ -1852,7 +2127,7 @@ The analytic solution for a general number of springs is given in Appendix \ref{ START_TEST(test_double_hooke) { double k1=5, k2=4, x=3, last_x=2.0, F, Fe, x1e, x2e; - one_dim_func_t *handlers[] = {&hooke, &hooke, NULL}; + one_dim_fn_t *handlers[] = {&hooke, &hooke, NULL}; void *data[] = {&k1, &k2, NULL}; double xi[] = {0.0, 0.0, 0.0}; int active[] = {1, 1, 0}; @@ -1896,6 +2171,12 @@ The interface is defined in [[list.h]], the implementation in [[list.c]], and th <> <> <> +<> +<> +<> +<> +<> +<> @ <>= @@ -1904,6 +2185,12 @@ The interface is defined in [[list.h]], the implementation in [[list.c]], and th <> <> <> +<> +<> +<> +<> +<> +<> @ <>= @@ -1925,8 +2212,9 @@ typedef struct list_struct { struct list_struct *prev; void *d; /* data */ } list_t; -@ +<> +@ [[head]] and [[tail]] return pointers to the head and tail nodes of the list: <>= @@ -1974,15 +2262,17 @@ int list_length(list_t *list) } @ -[[push]] inserts a node after the current position in the list +[[push]] inserts a node relative to the current position in the list without changing the current position. However, if the list we're pushing onto is [[NULL]], the current position isn't defined so we set it to that of the pushed domain. +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); +void push(list_t **pList, void *data, int below); @ <>= -void push(list_t **pList, void *data) +void push(list_t **pList, void *data, int below) { list_t *list, *new_node; assert(pList != NULL); @@ -1990,7 +2280,13 @@ void push(list_t **pList, void *data) new_node = create_node(data); if (list == NULL) *pList = new_node; - else { + else if (below==0) { /* insert above */ + if (list->prev != NULL) + list->prev->next = new_node; + new_node->prev = list->prev; + list->prev = new_node; + new_node->next = list; + } else { /* insert below */ if (list->next != NULL) list->next->prev = new_node; new_node->next = list->next; @@ -2000,10 +2296,9 @@ void push(list_t **pList, void *data) } @ -[[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. +[[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); @ @@ -2030,6 +2325,175 @@ void *pop(list_t **pList) } @ +[[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) +{ + list_t *list; + void *data; + assert(pList != NULL); + list = *pList; + *pList = NULL; + assert(list != NULL); /* not an empty list */ + while (list != NULL) { + data = pop(&list); + if (destroy != NULL) + 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) +{ + void **array; + int i,length; + assert(list != NULL); + assert(pLength != NULL); + assert(pArray != NULL); + + length = list_length(list); + array = (void **)malloc(sizeof(void *)*length); + assert(array != NULL); + list = head(list); + i=0; + while (list != NULL) { + array[i++] = list->d; + list = list->next; + } + *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) +{ + list_t *list, *flipper; + void *data; + assert(pList != NULL); + list = *pList; /* a>B>c>d */ + assert(list != NULL); /* not an empty list */ + if (downstream == 0) + flipper = list->prev; /* flipper is A */ + else + flipper = list->next; /* flipper is C */ + /* ensure there is somewhere to go in the direction we're heading */ + assert(flipper != NULL); + /* remove the current element */ + data = pop(&list); /* data=B, a>C>d */ + /* and put it back in in the correct spot */ + list = flipper; + if (downstream == 0) { + push(&list, data, 0); /* 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) +{ + int a,b; + a = *((int *)A); + b = *((int *)B); + if (a > b) return 1; + 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) +{ + list_t *list; + int stable=0; + assert(pList != NULL); + list = *pList; + assert(list != NULL); /* not an empty list */ + while (stable == 0) { + stable = 1; + list = head(list); + while (list->next != NULL) { + if ((*comp)(list->d, list->next->d) > 0) { + stable = 0; + move(&list, 1 /* downstream */); + } else + list = list->next; + } + } + *pList = list; +} +@ + +[[unique]] removes duplicates from a sorted list according to a user +specified comparison. Don't do this unless you have other pointers +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) +{ + list_t *list; + assert(pList != NULL); + list = *pList; + assert(list != NULL); /* not an empty list */ + list = head(list); + while (list->next != NULL) { + if ((*comp)(list->d, list->next->d) == 0) { + pop(&list); + } else + list = list->next; + } + *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) +{ + fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name); + list = head(list); + while (list != NULL) { + fprintf(file, " %p", list->d); + list = list->next; + } + fprintf(file, "\n"); +} +@ + [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]]. <>= list_t *create_node(void *data) @@ -2062,7 +2526,8 @@ The user must free the data pointed to by the node on their own. <>= #include /* assert() */ #include /* malloc(), free(), rand() */ -//#include /* fprintf(), stdout */ +#include /* fprintf(), stdout, FILE */ +#include "global.h" /* destroy_data_func_t */ @ \subsection{List unit tests} @@ -2105,7 +2570,7 @@ START_TEST(test_push) list_t *list=NULL; int i, p, e, n=3; for (i=0; id == 0 ); for (i=0; i>= <> <> -<> <> <> <> <> +<> +<> @ <>= @@ -2622,7 +3088,7 @@ clean_tension_model_utils : For unstretchable domains. <>= -{TG_NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL} +{NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL} @ \subsection{Constant} @@ -2713,7 +3179,7 @@ char const_tension_param_string[] = "0"; @ <>= -{TG_CONST, "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, "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} @@ -2807,7 +3273,7 @@ char hooke_param_string[]="0.05"; @ <>= -{TG_HOOKE, "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, "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} @@ -2935,7 +3401,7 @@ char wlc_param_string[]="0.39e-9,28.4e-9"; @ <>= -{TG_WLC, "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, "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. @@ -3057,7 +3523,7 @@ char fjc_param_string[]="0.5e-9,200"; @ <>= -{TG_FJC, "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, "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} @@ -3089,6 +3555,7 @@ char fjc_param_string[]="0.5e-9,200"; @ <>= +<> <> <> <> @@ -3130,18 +3597,18 @@ int main(int argc, char **argv) void *params; environment_t env; unsigned int flags; - one_dim_func_t *tension_handler[NUM_TENSION_GROUPS] = {0}; + one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0}; tension_handler_data_t tdata; int num_param_args; /* for INIT_MODEL() */ char **param_args; /* for INIT_MODEL() */ - get_options(argc, argv, &env, NUM_TENSION_GROUPS, tension_models, &model, &flags); + get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags); setup(tension_handler); if (flags & VFLAG) { printf("#initializing model %s with parameters %s\n", model->name, model->params); } INIT_MODEL("utility", model, params); tdata.group = NULL; - push(&tdata.group, params); + push(&tdata.group, params, 1); tdata.env = &env; tdata.persist = NULL; if (tension_handler[model->tg_group] == NULL) { @@ -3234,7 +3701,7 @@ void get_options(int argc, char **argv, environment_t *env, { char *prog_name = NULL; char c, options[] = "T:C:m:a:Vh"; - int tension_model=0; + int tension_model=0, num_strings; extern char *optarg; extern int optind, optopt, opterr; @@ -3252,6 +3719,7 @@ void get_options(int argc, char **argv, environment_t *env, case 'T': env->T = atof(optarg); break; case 'C': env->T = atof(optarg)+273.15; break; case 'm': + parse_list_string(optarg, ',', NULL, NULL, &num_strings, string_array); tension_model = index_tension_model(n_tension_models, tension_models, optarg); *model = tension_models+tension_model; break; @@ -4268,7 +4736,7 @@ END_TEST START_TEST(test_something) { double k=5, x=3, last_x=2.0, F; - one_dim_func_t *handlers[] = {&hooke}; + one_dim_fn_t *handlers[] = {&hooke}; void *data[] = {&k}; double xi[] = {0.0}; int active[] = {1}; -- 2.26.2