\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}
\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}
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.
<<version definition>>=
#define VERSION "0.4"
\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}),
\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}),
\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.
<<main program>>=
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;
}
[[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.
<<environment definition>>=
typedef struct environment_struct {
double T;
\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.
+
<<find tension definitions>>=
-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
+
+@
+
+<<tension handler array global declaration>>=
+extern one_dim_fn_t *tension_handlers[];
+@
+
+<<tension handler array global>>=
+one_dim_fn_t *tension_handlers[] = {
+ NULL,
+ &const_tension_handler,
+ &hooke_handler,
+ &wlc_handler,
+ &fjc_handler,
+ };
+
+@
+<<tension model global declarations>>=
+<<tension handler array global declaration>>
@
<<tension handler types>>=
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]].
<<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; i<NUM_TENSION_GROUPS; i++) {
- pdata[i] = data+i;
- data[i].group = get_group_list(domain_list, i);
- data[i].env = env;
- /* data[i].persist continues unchanged */
- if (data[i].group && i != TG_NULL) {
- /* an active group */
- if (active_groups[i] == 0) /* newly active */
- new_active_groups = 1;
- active_groups[i] = 1;
- } else {
- if (active_groups[i] == 1) /* newly inactive */
- new_active_groups = 1;
- active_groups[i] = 0;
+#ifdef DEBUG
+ fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
+ list_print(stderr, domain_list, "domain list");
+#endif
+
+ if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
+ /* free old statics */
+ if (per_group_handlers != NULL)
+ free(per_group_handlers);
+ if (per_group_params != NULL) {
+ for (i=0; i < num_active_groups; i++) {
+ assert(per_group_params[i] != NULL);
+ free(per_group_params[i]);
+ }
+ free(per_group_params);
}
- }
- F = tension_balance(NUM_TENSION_GROUPS, tension_handler, pdata, active_groups, new_active_groups, xi, last_x, x);
+
+ /* get new statics */
+ get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
+
+ /* fill in the group handlers and params */
+ for (i=0; i<num_active_groups; i++) {
+ tdata = (tension_handler_data_t *) per_group_params[i];
+#ifdef DEBUG
+ fprintf(stderr, "active group %d of %d", i, num_active_groups);
+ list_print(stderr, tdata->group, " 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}
$$
\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.
$$
}
@ [[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.
<<randomly unfold domains>>=
-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))) {
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}.
\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.
<<get options definitions>>=
<<model getopt definitions>>
\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.
<<tension model getopt definitions>>=
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 array definition>>=
-tension_model_getopt_t tension_models[NUM_TENSION_GROUPS] = {
+tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
<<null tension model getopt>>,
<<constant tension model getopt>>,
<<hooke tension model getopt>>,
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);
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);
}
@
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; i<n_tension_models; i++) {
+ fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
+ i, tension_models[i].name, tension_models[i].handler);
+ assert(tension_models[i].handler == tension_handlers[i]);
+ }
+#endif
+
/* setup defaults */
flags = FLAG_OUTPUT_UNFOLDING_FORCES;
prog_name = argv[0];
case 'T': env->T = 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);
n = atoi(optarg);
assert(n > 0);
for (i=0; i<n; i++) {
- push(pDomain_list, generate_domain(DS_FOLDED, tension_models+ftension_model, tension_models+utension_model, k_models+k_model));
+ push(pDomain_list, generate_domain(DS_FOLDED,
+ tension_models+ftension_model,
+ ftension_group,
+ ftension_params,
+ tension_models+utension_model,
+ utension_group,
+ utension_params,
+ k_models+k_model), 1);
}
*pNum_folded += n;
break;
n = atoi(optarg);
assert(n > 0);
for (i=0; i<n; i++) {
- push(pDomain_list, generate_domain(DS_UNFOLDED, tension_models+ftension_model, tension_models+utension_model, k_models+k_model));
+ push(pDomain_list, generate_domain(DS_UNFOLDED,
+ tension_models+ftension_model,
+ ftension_group,
+ ftension_params,
+ tension_models+utension_model,
+ utension_group,
+ utension_params,
+ k_models+k_model), 1);
}
break;
case 'V':
* INIT_MODEL("folded", folded_model, folded_params);
* defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
*/
-#define INIT_MODEL(role, model, param_pointer) \
- do { \
- parse_list_string(model->params, 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);
@
<<generate domain>>=
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
+ );
}
@
@
<<setup>>=
-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;
}
@
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.
<<search dt>>=
-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)
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);
/* 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);
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);
/* 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;
we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
<<determine dt>>=
<<search dt>>
-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.
//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;
\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.
<<domain definitions>>=
enum domain_state_t {DS_FOLDED,
DS_UNFOLDED
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 */
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)
{
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;
}
}
@
-\subsection{Group handling}
+\subsection{Domain model and group handling}
<<group functions>>=
+<<get tension handler>>
<<get group>>
+<<int comparison function>>
+<<find possible groups>>
+<<is group active>>
<<get group list>>
+<<get active groups>>
+@
+
+<<get tension handler>>=
+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;
+ }
+}
@
<<get group>>=
-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.
+<<find possible groups>>=
+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.
+<<is group 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.
<<get group list>>=
-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.
+<<get 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; i<num_tension_handlers; i++) {
+ handler = tension_handlers[i];
+ if (handler == NULL) continue; /* tensionless handler */
+ possible_group_numbers = head(find_possible_groups(list, handler));
+ while (possible_group_numbers != NULL) {
+ group = (int *)pop(&possible_group_numbers);
+ if (is_group_active(list, handler, *group) == 1) {
+ active_group_list = get_group_list(list, handler, *group);
+ push(&active_groups_list, active_group_list, 1);
+ push(&per_group_handlers_list, handler, 1);
+#ifdef DEBUG
+ fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
+ 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; i<num_active_groups; i++) {
+ handler = pop(&per_group_handlers_list);
+ per_group_handlers[i] = handler;
+
+ active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
+ assert(active_groups[i] != NULL);
+ active_group_list = pop(&active_groups_list);
+ ((tension_handler_data_t *)active_groups[i])->group = 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}
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.
<<parse delimited list string functions>>=
/* TODO, split out into parse.hc */
static int next_delim_index(char *string, char sep, char deeper, char shallower)
\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]]
<<tension-balance.h>>=
<<license comment>>
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.
<<tension balance function declaration>>=
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);
@
<<tension balance function>>=
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<num_tension_groups; i++) { /* remove null groups */
- if (active[i]) {
- active_handlers[active_groups] = tension_handler[i];
- active_params[active_groups] = params[i];
- active_xi[active_groups] = xi[i]; /* last balanced x_i if no new active groups */
- active_groups += 1;
- }
- }
- x_xo_data.n_groups = active_groups;
- x_xo_data.tension_handler = active_handlers;
- x_xo_data.handler_data = active_params;
- x_xo_data.xi = active_xi;
- if (active_groups == 1) { /* only one group, no balancing required */
+ assert(tension_handlers != NULL);
+ assert(params != NULL);
+ assert(num_tension_groups > 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<num_tension_groups; i++) printf(" %d: %g\t", i, xi[i]);
+ //for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, x_xo_data.xi[i]);
//printf("\n");
- if (new_active_groups || 1) { /* take a rough guess and attempt boundary conditions. */
- //printf("search bracket\n");
- oneD_bracket(x_of_xo, &x_xo_data, x, x/active_groups, &lb, &ub);
+ if (last_x == -1) { /* new group setup, reset x_xo_data */
+ /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
+ if (x_xo_data.xi != NULL)
+ free(x_xo_data.xi);
+ /* construct new x_xo_data */
+ x_xo_data.n_groups = num_tension_groups;
+ x_xo_data.tension_handler = tension_handlers;
+ x_xo_data.handler_data = params;
+ x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
+ for (i=0; i<num_tension_groups; i++)
+ x_xo_data.xi[i] = -1.0;
+ }
+ if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
+ if (x == 0) x_guess = 1.0;
+ else x_guess = x/num_tension_groups;
+#ifdef DEBUG
+ fprintf(stderr, "search bracket for %g with guess of %g\n", x, x_guess);
+#endif
+ oneD_bracket(x_of_xo, &x_xo_data, x, x_guess, &lb, &ub);
} else { /* work off the last balanced point */
- //printf("step off the old bracketing x %g + %g = %g\n", last_x, x-last_x, x);
+#ifdef DEBUG
+ fprintf(stderr, "step off the old bracketing x %g + %g = %g\n", last_x, x-last_x, x);
+#endif
if (x > 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<num_tension_groups; i++) {
- if (active[i]) {
- xi[i] = active_xi[active_groups];
- active_groups += 1;
+ if (num_tension_groups > 1 && 0) {
+ printf("balanced for x = %g with ", x);
+ for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, x_xo_data.xi[i]);
+ printf("\n");
}
}
-
- if (active_groups > 1 && 0) {
- printf("balanced for x = %g with ", x);
- for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, xi[i]);
- printf("\n");
- }
- free(active_handlers);
- free(active_params);
- free(active_xi);
-
+ F = (*tension_handlers[0])(xo, params[0]);
return F;
}
@
<<x of x0 definitions>>=
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;
@
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;
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;
Simple bisection, so it's robust and converges fairly quickly.
<<one dimensional function definition>>=
/* 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);
@
<<one dimensional solve declaration>>=
-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);
@
<<one dimensional solve>>=
/* 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)
{
The one dimensional solver needs a bracketed solution, which sometimes we don't have.
Generate bracketing $x$ values through bisection or exponential growth.
<<one dimensional bracket declaration>>=
-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);
@
<<one dimensional bracket>>=
-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;
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};
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};
<<list length declaration>>
<<push declaration>>
<<pop declaration>>
+<<list destroy declaration>>
+<<list to array declaration>>
+<<move declaration>>
+<<sort declaration>>
+<<unique declaration>>
+<<list print declaration>>
@
<<list functions>>=
<<list length>>
<<push>>
<<pop>>
+<<list destroy>>
+<<list to array>>
+<<move>>
+<<sort>>
+<<unique>>
+<<list print>>
@
<<list module makefile lines>>=
struct list_struct *prev;
void *d; /* data */
} list_t;
-@
+<<comparison function typedef>>
+@
[[head]] and [[tail]] return pointers to the head and tail nodes of the list:
<<head and tail declarations>>=
}
@
-[[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).
<<push declaration>>=
-void push(list_t **pList, void *data);
+void push(list_t **pList, void *data, int below);
@
<<push>>=
-void push(list_t **pList, void *data)
+void push(list_t **pList, void *data, int below)
{
list_t *list, *new_node;
assert(pList != NULL);
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;
}
@
-[[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.
<<pop declaration>>=
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.
+<<list destroy declaration>>=
+void list_destroy(list_t **pList, destroy_data_func_t *destroy);
+@
+<<list 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.
+<<list to array declaration>>=
+void list_to_array(list_t *list, int *length, void ***array);
+@
+<<list to 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.
+<<move declaration>>=
+void move(list_t **pList, int downstream);
+@
+<<move>>=
+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.
+<<comparison function typedef>>=
+typedef int (comparison_fn_t)(void *A, void *B);
+@
+For example
+<<int comparison function>>=
+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;
+}
+@
+<<sort declaration>>=
+void sort(list_t **pList, comparison_fn_t *comp);
+@
+Here we use the bubble sort algorithm.
+<<sort>>=
+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.
+<<unique declaration>>=
+void unique(list_t **pList, comparison_fn_t *comp);
+@
+<<unique>>=
+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.
+<<list print declaration>>=
+void list_print(FILE *file, list_t *list, char *list_name);
+@
+<<list print>>=
+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]].
<<create and destroy node>>=
list_t *create_node(void *data)
<<list includes>>=
#include <assert.h> /* assert() */
#include <stdlib.h> /* malloc(), free(), rand() */
-//#include <stdio.h> /* fprintf(), stdout */
+#include <stdio.h> /* fprintf(), stdout, FILE */
+#include "global.h" /* destroy_data_func_t */
@
\subsection{List unit tests}
list_t *list=NULL;
int i, p, e, n=3;
for (i=0; i<n; i++)
- push(&list, (void *)i);
+ push(&list, (void *)i, 1);
fail_unless(list == head(list), NULL);
fail_unless( (int)list->d == 0 );
for (i=0; i<n; i++) {
<<tension-model.h>>=
<<license comment>>
<<tension handler types>>
-<<find tension definitions>>
<<constant tension model declarations>>
<<hooke tension model declarations>>
<<worm-like chain tension model declarations>>
<<freely-jointed chain tension model declarations>>
+<<find tension definitions>>
+<<tension model global declarations>>
@
<<tension model module makefile lines>>=
For unstretchable domains.
<<null tension model getopt>>=
-{TG_NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
+{NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
@
\subsection{Constant}
@
<<constant tension model getopt>>=
-{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}
@
<<hooke tension model getopt>>=
-{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}
@
<<worm-like chain tension model getopt>>=
-{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.
@
<<freely-jointed chain tension model getopt>>=
-{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}
@
<<tension model globals>>=
+<<tension handler array global>>
<<constant tension model globals>>
<<hooke tension model globals>>
<<worm-like chain tension model globals>>
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) {
{
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;
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;
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};
projects / sawsim.git / commitdiff