Fix omega -> \omega in manual's timescale discussion.
[sawsim.git] / src / sawsim.nw
1 % sawsim - simulating single molecule mechanical unfolding.
2 % Copyright (C) 2008-2010  William Trevor King
3 %
4 % This program is free software: you can redistribute it and/or modify
5 % it under the terms of the GNU General Public License as published by
6 % the Free Software Foundation, either version 3 of the License, or
7 % (at your option) any later version.
8 %
9 % This program is distributed in the hope that it will be useful,
10 % but WITHOUT ANY WARRANTY; without even the implied warranty of
11 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 % GNU General Public License for more details.
13 %
14 % You should have received a copy of the GNU General Public License
15 % along with this program.  If not, see <http://www.gnu.org/licenses/>.
16 %
17 % The author may be contacted at <wking@drexel.edu> on the Internet, or
18 % write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
19 % Philadelphia PA 19104, USA.
20
21 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
22 % we have our own preamble,
23 % use `noweave -delay` to not print noweb's automatic one
24 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
25 \documentclass[letterpaper, 11pt]{article}
26 \usepackage{noweb}
27 \pagestyle{noweb}
28 \noweboptions{smallcode}
29
30 \usepackage{url}    % needed for natbib for underscores in urls?
31 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
32 % breaklinks breaks long links
33 % colorlinks colors link text instead of boxing it
34 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
35 \usepackage{amsthm}  % for theorem and proof environments
36 \newtheorem{theorem}{Theorem}
37
38 \usepackage{subfig} % Support side-by-side figures
39 \usepackage{pgf}    % Fancy graphics
40 \usepackage{tikz}   % A nice, inline PGF frontend
41 \usetikzlibrary{automata} % Graph-theory library
42
43 \usepackage{doc}    % BibTeX
44 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
45 % super selects citations in superscript mode
46 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
47
48 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
49 \usepackage[utf8x]{inputenc}  % allow UTF-8 characters
50
51 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
52 %\bibliographystyle{plain} % pick the bibliography style, includes dates
53 \bibliographystyle{plainnat}
54 \defcitealias{sw:noweb}{{\tt noweb}}
55 \defcitealias{galassi05}{GSL}
56 \defcitealias{sw:check}{{\tt check}}
57 \defcitealias{sw:python}{Python}
58
59 \topmargin -1.0in
60 \headheight 0.5in
61 \headsep 0.5in
62 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
63 \oddsidemargin -0.5in
64 \textwidth 7.5in
65 \setlength{\parindent}{0pt}
66 \setlength{\parskip}{5pt}
67
68 % For some reason, the \maketitle wants to go on it's own page
69 % so we'll just hardcode the following in our first page.
70 %\title{Sawsim: a sawtooth protein unfolding simulator}
71 %\author{W.~Trevor King}
72 %\date{\today}
73
74 \newcommand{\prog}{[[sawsim]]}
75 \newcommand{\lang}{[[C]]}
76 \newcommand{\p}[3]{\left#1 #2 \right#3}   % e.g. \p({complicated expr})
77 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}}   % ordinal th
78 \newcommand{\U}[1]{\ensuremath{\text{ #1}}}      % units: $1\U{m/s}$
79
80 % single spaced lists, from Alvin Alexander
81 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
82 \newenvironment{packed_item}{
83 \begin{itemize}
84   \setlength{\itemsep}{1pt}
85   \setlength{\parskip}{0pt}
86   \setlength{\parsep}{0pt}
87 }{\end{itemize}}
88
89 \begin{document}
90 \nwfilename{sawsim.nw}
91 \nwbegindocs{0}
92
93 %\maketitle
94 \begin{centering}
95   Sawsim: a sawtooth protein unfolding simulator \\
96   W.~Trevor King \\
97   \today \\
98 \end{centering}
99 \bigskip
100
101 \tableofcontents
102
103 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
104
105 \section{Introduction}
106
107 The unfolding of multi-globular protein strings is a stochastic
108 process.  In the \prog\ program, we use Monte Carlo methods to
109 simulate this unfolding through an explicit time-stepping approach.
110 We develop a program in \lang\ to simulate probability of unfolding
111 under a constant extension velocity of a chain of globular domains.
112
113 In order to extract physical meaning from many mechanical unfolding
114 simulations, it is necessary to compare the experimental results
115 against the results of simulations using different models and
116 parameters.  The model and parameters that best match the experimental
117 data provide the ‘best’ model for that experiment.
118
119 Previous mechanical unfolding studies have rolled their own
120 simulations for modelling their experiment, and there has been little
121 exchange and standardization of these simulations between groups.
122 The \prog\ program attempts to fill this gap, providing a flexible,
123 extensible, \emph{community} program, to make it easier for groups
124 to share the burden and increase the transparency of data analysis.
125
126 ‘Sawsim’ is blend of ‘sawtooth simulation’.
127
128 \subsection{Related work}
129
130 Sawim has been published\citep{king10}!  It is, as far as I know, the
131 only open-source Monte Carlo force spectroscopy simulator, but similar
132 closed source simulators have been around since the seminal paper by
133 \citet{rief97a}.  In chronological order, major contributions have
134 come from
135
136 \begin{packed_item}
137   \item \citet{rief97a} \citeyear{rief97a}:
138   \item \citet{rief97b} \citeyear{rief97b}:
139   \item \citet{kellermayer97} \citeyear{kellermayer97}:
140   \item \citet{rief98} \citeyear{rief98}:
141     first paper to focus mainly on the simulation
142   \item \citet{oberhauser98} \citeyear{oberhauser98}:
143   \item \citet{carrion-vazquez99a} \citeyear{carrion-vazquez99a}:
144   \item \citet{best02} \citeyear{best02}:
145     pointed out large fit valley
146   \item \citet{zinober02} \citeyear{zinober02}:
147     introduced the scaffold effect
148   \item \citet{jollymore09} \citeyear{jollymore09}:
149   \item \citet{king10} \citeyear{king10}:
150     introduced sawsim
151 \end{packed_item}
152
153
154 \subsection{About this document}
155
156 This document was written in \citetalias{sw:noweb} with code blocks in
157 \lang\ and comment blocks in \LaTeX.
158
159 \subsection{System Requirements and Dependencies}
160
161
162 \subsection{Change Log}
163
164 Version 0.0 used only the unfolded domain WLC to determine the tension
165 and had a fixed timestep $dt$.
166
167 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
168 probability for a given domain was constant.
169
170 Version 0.2 added Kramers' model unfolding rate calculations, and a
171 switch to select Bell's or Kramers' model.  This induced a major
172 rewrite, introducing the tension group abstraction, and a split to
173 multiple \lang\ source files.  Also added a unit-test suites for
174 sanity checking, and switched to SI units throughout.
175
176 Version 0.3 added integrated Kramers' model unfolding rate
177 calculations.  Also replaced some of my hand-coded routines with GNU
178 Scientific Library \citepalias{galassi05} calls.
179
180 Version 0.4 added Python evaluation for the saddle-point Kramers'
181 model unfolding rate calculations, so that the model functions could
182 be controlled from the command line.  Also added interpolating lookup
183 tables to accelerate slow unfolding rate model evaluation, and command
184 line control of tension groups.
185
186 Version 0.5 added the stiffness environmental parameter and it's
187 associated unfolding models.
188
189 Version 0.6 generalized from two state unfolding to arbitrary
190 $N$-state rate reactions.  This allows simulations of
191 unfolding-refolding, metastable intermediates, etc., but required
192 reorganizing the command line interface.
193
194 Version 0.7 added tension model inverses, which seems to reduce
195 computation time by about a factor of 2.  I was expecting more, but
196 I'll take what I can get.
197
198 Version 0.8 fixed a minor bug in unfolding probability for
199 multi-domain groups.  The probability of at least one domain unfolding
200 had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
201 However, for small $P$ the two are equivalent.
202
203 Version 0.9 added the [[-P]] option to sawsim, setting the target
204 probability for the per-state transition $P_N$, not the per-domain
205 transisiton $P_1$.
206
207 Version 0.10 added the [[-F]] option to sawsim, setting a limit on the
208 force change $dF$ in a single timestep for continuous pulling
209 simulations.  It also added the piston tension model (Section
210 \ref{sec.piston_tension}), and adjusted the stiffness calculations to
211 handle domains with undefined stiffness.
212
213 <<version definition>>=
214 #define VERSION "0.10"
215 @
216
217 \subsection{License}
218
219 <<license>>=
220  sawsim - simulating single molecule mechanical unfolding.
221  Copyright (C) 2008-2010, William Trevor King
222
223  This program is free software; you can redistribute it and/or
224  modify it under the terms of the GNU General Public License as
225  published by the Free Software Foundation; either version 3 of the
226  License, or (at your option) any later version.
227
228  This program is distributed in the hope that it will be useful, but
229  WITHOUT ANY WARRANTY; without even the implied warranty of
230  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
231  See the GNU General Public License for more details.
232
233  You should have received a copy of the GNU General Public License
234  along with this program.  If not, see <http://www.gnu.org/licenses/>.
235
236  The author may be contacted at <wking@drexel.edu> on the Internet, or
237  write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
238  Philadelphia PA 19104, USA.
239 @
240
241 <<license comment>>=
242 /*
243  <<license>>
244  */
245
246 @
247
248 \section{Main}
249 \label{sec.main}
250
251 The unfolding system is stored as a chain of \emph{domains} (Figure
252 \ref{fig.domain_chain}).  Domains can be folded, globular protein
253 domains, unfolded protein linkers, AFM cantilevers, or other
254 stretchable system components.  The system is modeled as graph of
255 possible domain states with transitions between them (Figure
256 \ref{fig.domain_states}).
257
258 \begin{figure}
259 \begin{center}
260 \subfloat[][]{\label{fig.domain_chain}
261 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
262   \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
263   \node[state] (A)              {domain 1};
264   \node[state] (B) [below of=A] {domain 2};
265   \node[state] (C) [below of=B] {.~.~.};
266   \node[state] (D) [below of=C] {domain $N$};
267   \node        (S) [below of=D] {Surface};
268   \node        (E) [above of=A] {};
269
270   \path[-] (A) edge (B)
271            (B) edge node [right] {Tension} (C)
272            (C) edge (D)
273            (D) edge (S);
274   \path[->,green] (A) edge node [right,black] {Extension} (E);
275 \end{tikzpicture}
276 }
277 \qquad
278 \subfloat[][]{\label{fig.domain_states}
279 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
280   \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
281   \node[state] (A)              {cantilever};
282   \node[state] (C) [below of=A] {transition};
283   \node[state] (B) [left of=C]  {folded};
284   \node[state] (D) [right of=C] {unfolded};
285
286   \path (B) edge [bend left] node [above] {$k_1$} (C)
287         (C) edge [bend left] node [below] {$k_1'$} (B)
288             edge [bend left] node [above] {$k_2$} (D)
289         (D) edge [bend left] node [below] {$k_2'$} (C);
290 \end{tikzpicture}
291 }
292 \caption{\subref{fig.domain_chain} Extending a chain of domains.  One end
293   of the chain is fixed, while the other is extended at a constant
294   speed.  The domains are coupled with rigid linkers, so the domains
295   themselves must stretch to accomodate the extension.
296   \subref{fig.domain_states} Each domain exists in a discrete state.  At
297   each timestep, it may transition into another state following a
298   user-defined state matrix such as this one, showing a metastable
299   transition state and an explicit ``cantilever'' domain.}
300 \end{center}
301 \end{figure}
302
303 Each domain contributes to the total tension, which depends on the
304 chain extension.  The particular model for the tension as a function
305 of extension is handled generally with function pointers.  So far the
306 following models have been implemented:
307 \begin{packed_item}
308  \item  Null (Appendix \ref{sec.null_tension}),
309  \item  Constant (Appendix \ref{sec.const_tension}),
310  \item  Hookean spring (Appendix \ref{sec.hooke}),
311  \item  Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
312  \item  Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
313  \item  Piston (Appendix \ref{sec.piston_tension}),
314 \end{packed_item}
315
316 The tension model and parameters used for a particular domain depend
317 on the domain's current state.  The transition rates between states
318 are also handled generally with function pointers, with
319 implementations for
320 \begin{packed_item}
321  \item  Null (Appendix \ref{sec.null_k}),
322  \item  Bell's model (Appendix \ref{sec.bell}),
323  \item  Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
324  \item  Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
325  \item  Kramers' saddle point model (Appendix \ref{sec.kramers}).
326 \end{packed_item}
327
328 The unfolding of the chain domains is modeled in two stages.  First
329 the tension of the chain is calculated.  Then the tension (assumed to
330 be constant over the length of the chain, see Section
331 \ref{sec.timescales}), is applied to each domain, and used to
332 calculate the probability of that domain changing state in the next
333 timestep $dt$.  Then the time is stepped forward, and the process is
334 repeated.  The simulation is complete when a pre-set time limit
335 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
336 <<main program>>=
337 int main(int argc, char **argv)
338 {
339   <<tension model getopt array definition>>
340   <<k model getopt array definition>>
341
342   /* variables initialized in get_options() */
343   list_t *state_list=NULL, *transition_list=NULL;
344   environment_t env={0};
345   double P, dF, t_max, dt_max, v, x_step;
346   state_t *stop_state=NULL;
347
348   /* variables used in the simulation loop */
349   double t=0, x=0, dt, F; /* time, distance, timestep, force */
350   int transition=1;  /* boolean marking a transition for tension and timestep optimization */
351
352   get_options(argc, argv, &P, &dF, &t_max, &dt_max, &v, &x_step,
353               &stop_state, &env, NUM_TENSION_MODELS, tension_models,
354               NUM_K_MODELS, k_models, &state_list, &transition_list);
355   setup();
356
357   if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
358   if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
359   while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
360     F = find_tension(state_list, &env, x, transition, 0);
361     if (x_step == 0)
362       dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, v, transition);
363     else
364       dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, 0, transition);
365     transition=random_transitions(transition_list, F, dt, &env);
366     if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
367     t += dt;
368     if (x_step == 0) {
369       x += v*dt;
370     } else {
371       if (dt == dt_max) { /* step completed */
372         x += x_step;
373         dt_max = x_step / v;
374       } else { /* still working on this step */
375         dt_max -= dt;
376       }
377     }
378   }
379   destroy_state_list(state_list);
380   destroy_transition_list(transition_list);
381   free_accels();
382   return 0;
383 }
384 @ The meat of the simulation is bundled into the three functions
385 [[find_tension]], [[determine_dt]], and [[random_transitions]].
386 [[find_tension]] is discussed in Section \ref{sec.find_tension},
387 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
388 [[random_transitions]] in Section \ref{sec.transition_rate}.
389
390 The stretched distance is extended in one of two modes depending on
391 [[xstep]].  If [[xstep == 0]], then the pulling is continuous, at
392 least within the limits of the inherent discretization of a time
393 stepped approach.  At any rate, the timestep is limited by the maximum
394 allowed timestep [[dt_max]] and the maximum allowed unfolding
395 probability in a given step [[P]].  In the second mode [[xstep > 0]],
396 and the end to end distance increases in discrete steps of that
397 length.  The time between steps is chosen to maintain an average
398 unfolding velocity [[v]].  This approach is less work to simulate than
399 smooth pulling and also closer to the reality of many experiments, but
400 it is practically impossible to treat analytically.  With both pulling
401 styles implemented, the effects of the discretization can be easily
402 calculated, bridging the gap between theory and experiment.
403
404 Environmental parameters important in determining reaction rates and
405 tensions (e.g.\ temperature, pH) are stored in a single structure to
406 facilitate extension to more complicated models in the future.
407 <<environment definition>>=
408 typedef struct environment_struct {
409   double T;
410   double stiffness;
411 } environment_t;
412 @
413
414 <<global.h>>=
415 #ifndef GLOBAL_H
416 #define GLOBAL_H
417
418 #define DOUBLE_PRECISION 1e-12
419
420 <<environment definition>>
421 <<one dimensional function definition>>
422 <<create/destroy definitions>>
423 <<model globals>>
424 #endif /* GLOBAL_H */
425 @
426 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
427 included multiple times.
428
429 \section{Simulation functions}
430
431 <<simulation functions>>=
432 <<find tension>>
433 <<determine dt>>
434 <<happens>>
435 <<does domain transition>>
436 <<randomly transition domains>>
437 @
438
439 \subsection{Tension}
440 \label{sec.find_tension}
441
442 Because the stretched system may be made up of several parts (folded
443 domains, unfolded domains, spring-like cantilever, \ldots), we will
444 assign the domains to states with tension models and groups.  The
445 models are used to determine the tension of all the domain in a given
446 state following some model (Hookean spring, WLC, \ldots).  The domains
447 are assumed to be commutative, so ordering is ignored.  The
448 interactions between the groups are assumed to be linear, but the
449 interactions between domains of the same group need not be.  Each
450 model has a tension handler function, which gives the tension $F$
451 needed to stretch a given group of domains a distance $x$.
452
453 Groups represent collections of states which the model should treat as
454 a single entity.  To understand the purpose of groups, consider a
455 system with two unfolded domain states (e.g.\ for two different
456 proteins) that were both modeled as WLCs.  If both unfolded states had
457 the same persistence length, it would be wise to assign them to the
458 same group.  This would reduce both the computational cost of
459 balancing the tension and the error associated with the inter-group
460 interaction linearity.  Note that group numbers only matter
461 \emph{within} model classes, since grouping domains with seperate
462 models doesn't make sense.  Within designated groups, the tension
463 parameters for each domain are still checked for compatibility, so if
464 you accidentally assign incompatible domains to the same group you
465 should get an appropriate error message.
466
467 <<find tension definitions>>=
468 #define NUM_TENSION_MODELS 6
469
470 @
471
472 <<tension handler array global declaration>>=
473 extern one_dim_fn_t *tension_handlers[];
474 @
475
476 <<tension handler array global>>=
477 one_dim_fn_t *tension_handlers[] = {
478               NULL,
479               &const_tension_handler,
480               &hooke_handler,
481               &wlc_handler,
482               &fjc_handler,
483               };
484
485 @
486
487 <<tension model global declarations>>=
488 <<tension handler array global declaration>>
489 @
490
491 <<tension handler types>>=
492 typedef struct tension_handler_data_struct {
493   list_t *group_tension_params; /* one item for each domain in the group */
494   environment_t *env;
495   void          *persist;
496 } tension_handler_data_t;
497
498 @
499
500 After sorting the chain into separate groups $G_i$, with tension
501 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
502 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
503 \forall i,j$.  Note that there may be several states within each
504 group.  [[find_tension]] is basically a wrapper to feed proper input
505 into the [[tension_balance]] function.  [[transition]] should be set
506 to 0 if there were no transitions since the previous call to
507 [[find_tension]] to support some optimizations that assume a small
508 increase in tension between steps (only true if no transition has
509 occured).  While we're messing with the tension handlers and if
510 [[const_env==0]], we also grab a measure of the current linker
511 stiffness for the environmental variable, which is needed by some of
512 the unfolding rate models (e.g.\ the linker-stiffness corrected Bell
513 model, Appendix \ref{sec.kbell}).
514 <<find tension>>=
515 double find_tension(list_t *state_list, environment_t *env, double x,
516                     int transition, int const_env)
517 { // TODO: !cont_env -> step_call, only alter env or statics if step_call==1
518   static int num_active_groups=0;
519   static one_dim_fn_t **per_group_handlers = NULL;
520   static one_dim_fn_t **per_group_inverse_handlers = NULL;
521   static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
522   static double last_x;
523   tension_handler_data_t *tdata;
524   double F, *pStiffness;
525   int i;
526
527 #ifdef DEBUG
528   fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
529 #endif
530
531   if (transition != 0 || num_active_groups == 0) { /* setup statics */
532     /* get new statics, freeing the old ones if they aren't NULL */
533     get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
534
535     /* fill in the group handlers and params */
536 #ifdef DEBUG
537     fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p, (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
538 #endif
539     for (i=0; i<num_active_groups; i++) {
540       tdata = (tension_handler_data_t *) per_group_data[i];
541 #ifdef DEBUG
542       fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
543       list_print(stderr, tdata->group_tension_params, " parameters");
544 #endif
545       tdata->env = env;
546       /* tdata->persist continues unchanged */
547     }
548     last_x = -1.0;
549   } /* else continue with the current statics */
550
551   if (const_env == 1)
552     pStiffness = NULL;
553   else
554     pStiffness = &env->stiffness;
555
556   F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, pStiffness);
557
558   last_x = x;
559
560   return F;
561 }
562 @ For the moment, we will restrict our group boundaries to a single
563 dimension, so $\sum_i x_i = x$, our total extension (it is this
564 restriction that causes the groups to interact linearly).  We'll also
565 restrict our extensions to all be positive.  With these restrictions,
566 the problem of balancing the tensions reduces to solving a system of
567 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
568 the number of active groups.  In general this can be fairly
569 complicated, but without much loss of practicality we can restrict
570 ourselves to strictly monotonically increasing, non-negative tension
571 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
572 simpler.  For example, it guarantees the existence of a unique, real
573 solution for finite forces.  See Appendix \ref{sec.tension_balance}
574 for the tension balancing implementation.
575
576 Each group must define a parameter structure if the tension model is
577 parametrized,
578  a function to create the parameter structure,
579  a function to destroy the parameter structure, and
580  a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
581 The tension-specific parameter structures are contained in the domain
582 group field.  For optimization, a group may also define an inverse
583 tension function as an optimization.  See the Section
584 \ref{sec.model_selection} for a discussion on the command line
585 framework and inverse function discussion.  See the worm-like chain
586 implementation for an example (Section \ref{sec.wlc}).
587
588 The major limitation of our approach is that all of our tension models
589 are Markovian (memory-less), which is not really a problem for the
590 chains (see \citet{evans99} for WLC thermalization rate calculations).
591
592 \subsection{Transition rate}
593 \label{sec.transition_rate}
594
595 Each state transition is modeled as unimolecular, first order reaction
596 $$
597   \text{State 1} \xrightarrow{k} \text{State 2}
598 $$
599 With the rate constant $k$ defined by
600 $$
601   \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
602 $$
603 So the probability of a given domain transitioning in a short time
604 $dt$ is given by
605 $$
606   P_1 = k \cdot dt.
607 $$
608
609 For a group of $N$ identical domains, and therefore identical
610 unfolding probabilities $P_1$, the probability of $n$ domains
611 unfolding in a given timestep is given by the binomial distribution
612 $$
613   P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^{(N-n)}.
614 $$
615 The probability that \emph{none} of these domains unfold is then
616 $$
617   P(0) = (1-P_1)^N,
618 $$
619 so the probability that \emph{at least one} domain unfolds is
620 $$
621   P_N = 1-(1-P_1)^N.
622 $$
623 Note that for small $P$,
624 $$
625   p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
626     \approx NP_1.
627 $$
628 This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
629 <<PN definition>>=
630 /* find multi-domain transition rate from the single domain rate */
631 #define PN(P,N) (1.0-pow(1.0-(P), (N)))
632
633 @
634
635 We can also define a macro to find the approprate $dt$ to achieve a
636 target $P_N$ with a given $k$ and $N$.
637 \begin{align}
638   P_N &= 1-(1-P_1)^N \\
639   1-P_1 &= (1-P_N)^{1/N} \\
640   P_1 &= 1-(1-P_N)^{1/N}
641 \end{align}
642 <<P1 definition>>=
643 #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
644 @
645
646 We take some time to discuss the meaning of $p(n>1)$
647 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
648
649 <<does domain transition>>=
650 int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
651 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
652   double k;
653   k = accel_k(transition->k, F, env, transition->k_params);
654   //(*transition->k)(F, env, domain->k_params);
655   //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
656   return happens(PN(k*dt, N_INIT(transition)));
657 }
658 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
659
660 <<randomly transition domains>>=
661 int random_transitions(list_t *transition_list, double tension,
662                        double dt, environment_t *env)
663 { /* returns 1 if there was a transition and 0 otherwise */
664   while (transition_list != NULL) {
665     if (T(transition_list)->initial_state->num_domains > 0
666         && domain_transitions(tension, dt, env, T(transition_list))) {
667       if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
668         fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
669       T(transition_list)->initial_state->num_domains--;
670       T(transition_list)->final_state->num_domains++;
671       return 1; /* our one domain has transitioned, stop looking for others */
672     }
673     transition_list = transition_list->next;
674   }
675   return 0;
676 }
677 @
678
679 \subsection{Timescales}
680 \label{sec.timescales}
681
682 The simulation assumes that chain equilibration is instantaneous
683 relative to the simulation timestep $dt$, so that tension is uniform
684 along the chain.  The quality of this assumption depends on your
685 particular chain.  For example, a damped spring thermalizes on a
686 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
687 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
688 frequency $\omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
689 and $b$ sets the drag $F_b = -bv$.  For our cantilevers $\tau$ is
690 on the order of milliseconds, which is longer than a timestep.
691 However, the approximation is still reasonable, because there is
692 rarely unfolding during the cantilever return.  You could, of course,
693 take cantilever, WLC, etc.\ response times into effect, but that
694 would significantly complicate a simulation that seems to work
695 well enough without bothering :p.  For WLC and FJC relaxation times,
696 see the Appendix of \citet{evans99} and \citet{kroy07}.
697
698 Besides assuming our timestep is much greater than equilibration
699 times, we also force it to be short enough so that only one domain may
700 unfold in a given timestep.  For an unfolding timescale of $dt_u =
701 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
702 of two domains unfolding in the same timestep is small (see
703 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
704 [[main()]] in Section \ref{sec.main} set by [[-P]] in
705 [[get_options()]] in Section \ref{sec.get_options}). This approach
706 breaks down as the adaptive timestep scheme approaches the transition
707 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
708 \citep{klimov00}, so this shouldn't be a problem. To reassure
709 yourself, you can ask the simulator to print the smallest timestep
710 that was used with TODO.
711
712 Even if the likelyhood of two domains unfolding in the same timestep
713 is small, if you run long enough simulations it will eventually occur.
714 In this case, we assume that $dt_t \ll dt$, so even if two domains
715 would unfold if we stuck to our unfolding rate $k$ for an entire
716 timestep $dt$, in ``reality'' only one of those domains unfolds.
717 Because we do not know when in the timestep the unfolding took place,
718 we move on to the next timestep without checking for further
719 unfoldings.  This ``unchecked timestep portion'' should not contribute
720 significant errors to the calculation, because if $dt$ was small
721 enough that unfolding was unlikely at the pre-unfolding force, it
722 should be even more unlikely at the post-unfolding force, especially
723 over only a fraction of the timestep.  The only dangerous cases for
724 the current implementation are those where unfolding intermediates are
725 much less stable than their precursor states, in which case an
726 unfolding event during the remaining timestep may actually be likely.
727 A simple workaround for such cases is to pick the value for [[P]] to
728 be small enough that the $dt \ll dt_u$ assumption survives even under
729 a transition populating the unstable intermediate.
730
731 \subsection{Adaptive timesteps}
732 \label{sec.adaptive_dt}
733
734 We'd like to pick $dt$ so the probability of changing state
735 (e.g.\ unfolding) in the next timestep is small.  If we don't adapt our
736 timestep, we also risk breaking our approximation that $F(x' \in [x,
737   x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
738 given timestep.  The simulation would have been accurate for
739 sufficiently small fixed timesteps, but adaptive timesteps will allow
740 us to move through low tension regions in fewer steps, leading to a
741 more efficient simulation.
742
743 Consider the two state folded $\rightarrow$ unfolded transition.
744 Because $F(x)$ is monotonically increasing between unfoldings,
745 excessively large timesteps will lead to erroneously small unfolding
746 rates (an thus, higher average unfolding force).
747
748 The actual adaptive timestep implementation is not particularly
749 interesting, since we are only required to reduce $dt$ to somewhere
750 below a set threshold, so I've removed it to Appendix
751 \ref{sec.adaptive_dt_implementation}.
752
753 \section{Models}
754
755 TODO: model intro.
756
757 The models provide the physics of the simulation, but the simulation
758 allows interchangeable models, and we are currently focusing on the
759 simulation itself, so we remove the actual model implementation to the
760 appendices.
761
762 The tension models are defined in Appendix \ref{sec.tension_models},
763 and unfolding models are defined in Appendix \ref{sec.k_models}.
764
765 <<model globals>>=
766 #define k_B   1.3806503e-23  /* J/K */
767
768 @
769
770
771 \section{Command line interface}
772
773 <<get option functions>>=
774 <<help>>
775 <<index tension model>>
776 <<index k model>>
777 <<get options>>
778 @
779
780 \subsection{Model selection}
781 \label{sec.model_selection}
782
783 The main difficulty with the command line interface in \prog\ is
784 developing an intuitive method for accessing the possibly large number
785 of available models.  We'll treat this generally by defining an array
786 of available models, containing their important parameters.
787
788 <<get options definitions>>=
789 <<model getopt definitions>>
790 @
791
792 <<create/destroy definitions>>=
793 typedef void *create_data_func_t(char **param_strings);
794 typedef void destroy_data_func_t(void *param_struct);
795 @
796
797 <<model getopt definitions>>=
798 <<tension model getopt definitions>>
799 <<k model getopt definitions>>
800 @
801
802 In order to choose models, we need a method of assembling a reaction
803 graph such as that in Figure \ref{fig.domain_states} and population
804 the graph with a set of domains.  First we declare the domain states
805 following
806 \begin{center}
807 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
808 \end{center}
809 e.g.
810 \begin{center}
811 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
812 \end{center}
813 This creates the state named [[unfolded]], using the WLC model and one
814 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
815 The tension parameters are then set to [[1e-8,4e-10]], which in the
816 case of the WLC are the contour and persistence lengths respectively.
817 Finally we populate the state by adding five domains to the state.
818 The name of the state is importent for identifying states later.
819
820 Once the domains have been created and populated, you can added
821 transition rates following
822 \begin{center}
823   [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
824 \end{center}
825 e.g.
826 \begin{center}
827   [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
828 \end{center}
829 This creates a reaction going from the [[folded]] state to the
830 [[unfolded]] state, with the rate constant given by the Bell model
831 (Appendix \ref{sec.bell}).  The unfolding parameters are then set to
832 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
833 unforced rate and transition state distance respectively.
834
835 \subsubsection{Tension}
836
837 To access the various tension models from the command line, we define
838 a structure that contains all the neccessary information about the
839 model.
840 <<tension model getopt definitions>>=
841 typedef struct tension_model_getopt_struct {
842   one_dim_fn_t        *handler;     /* fn implementing the model on a group */
843   one_dim_fn_t        *inverse_handler; /* fn implementing inverse on group */
844   char                *name;        /* name string identifying the model    */
845   char                *description; /* a brief description of the model     */
846   int                  num_params;  /* number of extra params for the model */
847   char               **param_descriptions; /* descriptions of extra params  */
848   char                *params;      /* default values for the extra params  */
849   create_data_func_t  *creator;     /* param-string -> model-data-struct fn */
850   destroy_data_func_t *destructor;  /* fn to free the model data structure  */
851 } tension_model_getopt_t;
852 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
853 bisection.  Obviously, this will be slower than computing an
854 analytical inverse and more prone to rounding errors, but it may be
855 easier if a closed-form inverse does not exist.
856
857 <<tension model getopt array definition>>=
858 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
859   <<null tension model getopt>>,
860   <<constant tension model getopt>>,
861   <<hooke tension model getopt>>,
862   <<worm-like chain tension model getopt>>,
863   <<freely-jointed chain tension model getopt>>,
864   <<piston tension model getopt>>
865 };
866 @
867
868 \subsubsection{Unfolding rate}
869
870 <<k model getopt definitions>>=
871 #define NUM_K_MODELS 6
872
873 typedef struct k_model_getopt_struct {
874   char *name;
875   char *description;
876   k_func_t *k;
877   int num_params;
878   char **param_descriptions;
879   char *params;
880   create_data_func_t *creator;
881   destroy_data_func_t *destructor;
882 } k_model_getopt_t;
883 @
884
885 <<k model getopt array definition>>=
886 k_model_getopt_t k_models[NUM_K_MODELS] = {
887   <<null k model getopt>>,
888   <<const k model getopt>>,
889   <<bell k model getopt>>,
890   <<kbell k model getopt>>,
891   <<kramers k model getopt>>,
892   <<kramers integ k model getopt>>
893 };
894 @
895
896 \subsection{help}
897
898 <<help>>=
899 void help(char *prog_name, double P, double dF,
900           double t_max, double dt_max, double v, double x_step,
901           state_t *stop_state, environment_t *env,
902           int n_tension_models, tension_model_getopt_t *tension_models,
903           int n_k_models, k_model_getopt_t *k_models,
904           int k_model, list_t *state_list)
905 {
906   state_t *state=NULL;
907   int i, j;
908
909   if (state_list != NULL)
910     state = S(tail(state_list));
911
912   printf("usage: %s [options]\n", prog_name);
913   printf("Version %s\n\n", VERSION);
914   printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
915
916   printf("Simulation options:\n");
917   printf("");
918   printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
919   printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
920   printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
921   printf("-P P\tMaximum transition probability for dt (currently %g)\n", P);
922   printf("-F dF\tMaximum change in force over dt.  Only relevant if dx > 0. (currently %g N)\n", dF);
923   printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
924   printf("-x dx\tPulling step size (currently %g m)\n", x_step);
925   printf("\tWhen dx = 0, the pulling is continuous.\n");
926   printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
927
928   printf("Environmental options:\n");
929   printf("-T T\tTemperature T (currently %g K)\n", env->T);
930   printf("-C T\tYou can also set the temperature T in Celsius\n");
931
932   printf("State creation options:\n");
933   printf("-s name,model[[,group],params]\tCreate a new state.\n");
934   printf("Once you've set up the state, you may populate it with domains\n");
935   printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
936
937   printf("Transition creation options:\n");
938   printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
939
940   printf("To double check your reaction graph, you can produce graphviz dot output\n");
941   printf("describing the current states and transitions.\n");
942   printf("-D\tPrint dot description of states and transitions and exit.\n");
943
944   printf("Simulation output mode options:\n");
945   printf("There are two output modes.  In standard mode, only the transition\n");
946   printf("events are printed.  For example:\n");
947   printf("  #Force (N)\tInitial state\tFinal state\n");
948   printf("  123.456e-12\tfolded\tunfolded\n");
949   printf("  ...\n");
950   printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
951   printf("  #Distance (m)\tForce (N)\tStiffness (N/m)\n");
952   printf("  0.001\t0.0023\t0.002\n");
953   printf("  ...\n");
954   printf("-V\tChange output to verbose mode.\n");
955   printf("-h\tPrint this help and exit.\n");
956   printf("\n");
957
958   printf("Tension models:\n");
959   for (i=0; i<n_tension_models; i++) {
960     printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
961     for (j=0; j < tension_models[i].num_params ; j++)
962       printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
963     printf("\t\tdefaults: %s\n", tension_models[i].params);
964   }
965   printf("Unfolding rate models:\n");
966   for (i=0; i<n_k_models; i++) {
967     printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
968     for (j=0; j < k_models[i].num_params ; j++)
969       printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
970     printf("\t\tdefaults: %s\n", k_models[i].params);
971   }
972
973   printf("Examples:\n");
974   printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded.  Stop once there are no more folded domains.\n");
975   printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s folded,null -N8 -s 'unfolded,wlc,{0.39e-9,28e-9}' -k 'folded,unfolded,bell,{3.3e-4,0.25e-9}' -q folded\n", prog_name);
976   printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded.  Stop after 0.25 seconds\n");
977   printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s 'folded,wlc,{1e-8,4e-10}' -N8 -s 'unfolded,wlc,1,{0.4e-9,28.1e-9}' -k 'folded,unfolded,bell,{5e-5,0.225e-9}' -t0.25\n", prog_name);
978 }
979 @
980
981 \subsection{Get options}
982 \label{sec.get_options}
983
984 <<get options>>=
985 <<safe strto*>>
986 void get_options(int argc, char **argv, double *pP, double *pDF,
987                  double *pT_max, double *pDt_max, double *pV,
988                  double *pX_step, state_t **pStop_state, environment_t *env,
989                  int n_tension_models, tension_model_getopt_t *tension_models,
990                  int n_k_models, k_model_getopt_t *k_models,
991                  list_t **pState_list, list_t **pTransition_list)
992 {
993   char *prog_name = NULL;
994   char c, options[] = "q:t:d:P:F:v:x:T:C:s:N:k:DVh";
995   int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
996   char *state_name=NULL;
997   int i, n, tension_group=0;
998   list_t *temp_list=NULL;
999   int num_strings;
1000   char **strings;
1001   void *model_params; /* for INIT_MODEL() */
1002   int num_param_args; /* for INIT_MODEL() */
1003   char **param_args;  /* for INIT_MODEL() */
1004   extern char *optarg;
1005   extern int optind, optopt, opterr;
1006
1007   assert(argc > 0);
1008
1009 #ifdef DEBUG
1010   for (i=0; i<n_tension_models; i++) {
1011     fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
1012             i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
1013     assert(tension_models[i].handler == tension_handlers[i]);
1014   }
1015 #endif
1016
1017   /* setup defaults */
1018   flags = FLAG_OUTPUT_TRANSITION_FORCES;
1019   prog_name = argv[0];
1020   *pStop_state = NULL;
1021   *pT_max = -1;     /* s, -1 removes this limit */
1022   *pDt_max = 0.001; /* s                        */
1023   *pP = 1e-3;       /* % pop per s              */
1024   *pDF = 1e-12;     /* N                        */
1025   *pV = 1e-6;       /* m/s                      */
1026   *pX_step = 0.0;   /* m                        */
1027   env->T = 300.0;   /* K                        */
1028   *pState_list = NULL;
1029   *pTransition_list = NULL;
1030
1031   while ((c=getopt(argc, argv, options)) != -1) {
1032     switch(c) {
1033     case 'q':
1034       if (STRMATCH(optarg, "-")) {
1035         *pStop_state = NULL;
1036       } else {
1037         temp_list = head(*pState_list);
1038         while (temp_list != NULL) {
1039           if (STRMATCH(S(temp_list)->name, optarg)) {
1040             *pStop_state = S(temp_list);
1041             break;
1042           }
1043           temp_list = temp_list->next;
1044         }
1045       }
1046       break;
1047     case 't':  *pT_max  = safe_strtod(optarg, "-t");         break;
1048     case 'd':  *pDt_max = safe_strtod(optarg, "-d");         break;
1049     case 'P':  *pP      = safe_strtod(optarg, "-P");         break;
1050     case 'F':  *pDF     = safe_strtod(optarg, "-F");         break;
1051     case 'v':  *pV      = safe_strtod(optarg, "-v");         break;
1052     case 'x':  *pX_step = safe_strtod(optarg, "-x");         break;
1053     case 'T':  env->T   = safe_strtod(optarg, "-T");         break;
1054     case 'C':  env->T   = safe_strtod(optarg, "-C")+273.15;  break;
1055     case 's':
1056       parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1057       assert(num_strings >= 2 && num_strings <= 4);
1058
1059       state_name = strings[0];
1060       if (state_index(*pState_list, state_name) != -1) {
1061         fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
1062         exit(1);
1063       }
1064       tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
1065       if (num_strings == 4)
1066         tension_group = safe_strtoi(strings[2], optarg);
1067       else
1068         tension_group = 0;
1069       if (num_strings >= 3)
1070         tension_models[tension_model].params = strings[num_strings-1];
1071 #ifdef DEBUG
1072       fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s, group %d\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group);
1073 #endif
1074       INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
1075       push(pState_list, create_state(state_name,
1076                                      tension_models[tension_model].handler,
1077                                      tension_models[tension_model].inverse_handler,
1078                                      tension_group,
1079                                      model_params,
1080                                      tension_models[tension_model].destructor,
1081                                      0), 1);
1082       *pState_list = tail(*pState_list);
1083       state_name = NULL;
1084       free_parsed_list(num_strings, strings);
1085       break;
1086     case 'N':
1087       n = safe_strtoi(optarg, "-N");
1088 #ifdef DEBUG
1089       fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1090 #endif
1091       S(*pState_list)->num_domains += n;
1092       break;
1093     case 'k':
1094       parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1095       assert(num_strings >= 3 && num_strings <= 4);
1096       initial_state = state_index(*pState_list, strings[0]);
1097       final_state = state_index(*pState_list, strings[1]);
1098       k_model = index_k_model(n_k_models, k_models, strings[2]);
1099 #ifdef DEBUG
1100       fprintf(stderr, "%s:\t%s -> %s reaction from state %s (%d) -> state %s (%d)\n", __FUNCTION__, optarg, strings[2], strings[0], initial_state, strings[1], final_state);
1101 #endif
1102       assert(initial_state != final_state);
1103       if (num_strings == 4)
1104         k_models[k_model].params = strings[3];
1105       INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1106       push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1107                                                list_element(*pState_list, final_state),
1108                                                k_models[k_model].k,
1109                                                model_params,
1110                                                k_models[k_model].destructor), 1);
1111       free_parsed_list(num_strings, strings);
1112       break;
1113     case 'D':
1114       print_state_graph(stdout, *pState_list, *pTransition_list);
1115       exit(0);
1116       break;
1117     case 'V':
1118       flags = FLAG_OUTPUT_FULL_CURVE;
1119       break;
1120     case '?':
1121       fprintf(stderr, "unrecognized option '%c'\n", optopt);
1122       /* fall through to default case */
1123     default:
1124       help(prog_name, *pP, *pDF, *pT_max, *pDt_max, *pV, *pX_step,
1125            *pStop_state, env, n_tension_models, tension_models,
1126            n_k_models, k_models, k_model, *pState_list);
1127       exit(1);
1128     }
1129   }
1130   /* check the arguments */
1131   assert(*pState_list != NULL); /* no domains! */
1132   assert(*pP > 0.0); assert(*pP < 1.0);
1133   assert(*pDt_max > 0.0);
1134   assert(*pV > 0.0);
1135   assert(env->T > 0.0);
1136
1137   crosslink_groups(*pState_list, *pTransition_list);
1138
1139   return;
1140 }
1141 @
1142
1143 <<index tension model>>=
1144 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1145 {
1146   int i;
1147   for (i=0; i<n_models; i++)
1148     if (STRMATCH(models[i].name, name))
1149       break;
1150   if (i == n_models) {
1151     fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1152     exit(1);
1153   }
1154   return i;
1155 }
1156 @
1157
1158 <<index k model>>=
1159 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1160 {
1161   int i;
1162   for (i=0; i<n_models; i++)
1163     if (STRMATCH(models[i].name, name))
1164       break;
1165   if (i == n_models) {
1166     fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1167     exit(1);
1168   }
1169   return i;
1170 }
1171 @
1172
1173 <<initialize model definition>>=
1174 /* requires int num_param_args and char **param_args in the current scope
1175  * usage:
1176  *  INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1177  * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1178  */
1179 #define INIT_MODEL(role, model, param_string, param_pointer) \
1180   do {                                                       \
1181     parse_list_string((param_string), SEP, '{', '}',         \
1182                       &num_param_args, &param_args);         \
1183     if (num_param_args != (model)->num_params) {             \
1184       fprintf(stderr,                                        \
1185               "%s model %s expected %d params,\n",           \
1186               (role), (model)->name, (model)->num_params);   \
1187       fprintf(stderr,                                        \
1188               "not the %d params in '%s'\n",                 \
1189               num_param_args, (param_string));               \
1190       assert(num_param_args == (model)->num_params);         \
1191     }                                                        \
1192     if ((model)->creator)                                    \
1193       param_pointer = (*(model)->creator)(param_args);       \
1194     else                                                     \
1195       param_pointer = NULL;                                  \
1196     free_parsed_list(num_param_args, param_args);            \
1197   } while (0);
1198 @
1199
1200 First we define some safe string parsing functions.  Currently these
1201 abort on any irregularity, but in the future they could print error
1202 messages, etc.  [[static]] because the functions are currently
1203 defined in each [[*.c]] file that uses them, but perhaps they should
1204 be moved to [[utils.h/c]] or some such instead\ldots
1205
1206 <<safe strto*>>=
1207 static int safe_strtoi(const char *s, const char *description) {
1208   char *endp = NULL;
1209   long int ret;
1210   assert(s != NULL);
1211   ret = strtol(s, &endp, 10);
1212   if (endp[0] != '\0') { //strlen(endp) > 0) {
1213     fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1214             endp, s, description);
1215     assert(1==0);  //strlen(endp) == 0);
1216   } if (endp == s) {
1217     fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1218             s, description);
1219     assert(endp != s);
1220   } else if (errno == ERANGE) {
1221     fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1222     assert(errno != ERANGE);
1223   }
1224   return (int) ret;
1225 }
1226
1227 static double safe_strtod(const char *s, const char *description) {
1228   char *endp = NULL;
1229   double ret;
1230   assert(s != NULL);
1231   ret = strtod(s, &endp);
1232   if (endp[0] != '\0') { //strlen(endp) > 0) {
1233     fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1234             endp, s, description);
1235     assert(1==0);  //strlen(endp) == 0);
1236   } if (endp == s) {
1237     fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1238             s, description);
1239     assert(endp != s);
1240   } else if (errno == ERANGE) {
1241     fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1242     assert(errno != ERANGE);
1243   }
1244   return ret;
1245 }
1246
1247 @
1248
1249 \phantomsection
1250 \appendix
1251 \addcontentsline{toc}{section}{Appendicies}
1252
1253 \section{sawsim.c details}
1254
1255 \subsection{Layout}
1256
1257 The general layout of our simulation code is:
1258 <<*>>=
1259 //#define DEBUG
1260 <<license comment>>
1261 <<includes>>
1262 <<definitions>>
1263 <<globals>>
1264 <<functions>>
1265 <<main program>>
1266 @
1267
1268 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1269 <<includes>>=
1270 #include <assert.h> /* assert()                                */
1271 #include <stdlib.h> /* malloc(), free(), rand(), strto*()      */
1272 #include <stdio.h>  /* fprintf(), stdout                       */
1273 #include <string.h> /* strlen, strtok()                        */
1274 #include <math.h>   /* exp(), M_PI, sqrt()                     */
1275 #include <sys/types.h> /* pid_t (returned by getpid())         */
1276 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1277 #include <errno.h>  /* errno, ERANGE (for safe_strto*())       */
1278 #include "global.h"
1279 #include "list.h"
1280 #include "tension_balance.h"
1281 #include "k_model.h"
1282 #include "tension_model.h"
1283 #include "parse.h"
1284 #include "accel_k.h"
1285
1286 @
1287
1288 <<definitions>>=
1289 <<version definition>>
1290 <<flag definitions>>
1291 <<max/min definitions>>
1292 <<string comparison definition and globals>>
1293 <<initialize model definition>>
1294 <<get options definitions>>
1295 <<state definitions>>
1296 <<transition definitions>>
1297 @
1298
1299 <<globals>>=
1300 <<flag globals>>
1301 <<model globals>>
1302 @
1303
1304 <<functions>>=
1305 <<create/destroy state>>
1306 <<destroy state list>>
1307 <<state index>>
1308 <<create/destroy transition>>
1309 <<destroy transition list>>
1310 <<print state graph>>
1311 <<group functions>>
1312 <<simulation functions>>
1313 <<boilerplate functions>>
1314 @
1315
1316 <<boilerplate functions>>=
1317 <<setup>>
1318 <<get option functions>>
1319 @
1320
1321 <<setup>>=
1322 void setup(void)
1323 {
1324   srand(getpid()*time(NULL)); /* seed rand() */
1325 }
1326 @
1327
1328 <<flag definitions>>=
1329 /* in octal b/c of prefixed '0' */
1330 #define FLAG_OUTPUT_FULL_CURVE        01
1331 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1332 @
1333
1334 <<flag globals>>=
1335 static unsigned long int flags = 0;
1336 @
1337
1338 \subsection{Utilities}
1339 \label{sec.utils}
1340
1341 <<max/min definitions>>=
1342 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1343 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1344 @
1345
1346 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1347 <<string comparison definition and globals>>=
1348 // Check if two strings match, return 1 if they do
1349 static char *temp_string_A;
1350 static char *temp_string_B;
1351 #define STRMATCH(a,b)   (temp_string_A=a, temp_string_B=b, \
1352         strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1353         !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1354         /* +1 to also compare the '\0' */
1355 @
1356
1357 We also define a macro for our [[check]] unit testing
1358 <<check relative error>>=
1359 #define CHECK_ERR(max_err, expected, received)                               \
1360   do {                                                                       \
1361     fail_unless( (received-expected)/expected < max_err,                     \
1362                  "relative error %g >= %g in %s (Expected %g, received %g)", \
1363                  (received-expected)/expected, max_err, #received,           \
1364                  expected, received);                                        \
1365     fail_unless(-(received-expected)/expected < max_err,                     \
1366                  "relative error %g >= %g in %s (Expected %g, received %g)", \
1367                 -(received-expected)/expected, max_err, #received,           \
1368                  expected, received);                                        \
1369   } while(0)
1370
1371 @
1372
1373 <<happens>>=
1374 int happens(double probability)
1375 {
1376   assert(probability >= 0.0); assert(probability <= 1.0);
1377   return (double)rand()/RAND_MAX < probability; /* TODO: replace with GSL rand http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html x*/
1378 }
1379 @
1380
1381
1382 \subsection{Adaptive timesteps}
1383 \label{sec.adaptive_dt_implementation}
1384
1385 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1386 chain model), so basing the timestep on the unfolding probability at
1387 the current tension is dangerous, and we need to search for a $dt$ for
1388 which $P_N(F(x+v*dt))<P_\text{max}$ (See Section
1389 \ref{sec.transition_rate} for a discussion of $P_N$).  There are two
1390 cases to consider.  In the most common, no domains have unfolded since
1391 the last step, and we expect the next step to be slightly shorter than
1392 the previous one.  In the less common, domains did unfold in the last
1393 step, and we expect the next step to be considerably longer than the
1394 previous one.
1395 <<search dt>>=
1396 double search_dt(transition_t *transition,
1397                  list_t *state_list,
1398                  environment_t *env, double x,
1399                  double max_prob, double max_dt, double v)
1400 { /* Returns a valid timestep dt in seconds for the given transition.
1401    * Takes a list of all states, a pointer env to the environmental data,
1402    * a starting separation x in m, a max_prob between 0 and 1,
1403    * max_dt in s, stretching velocity v in m/s.
1404    */
1405   double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
1406
1407   /* get upper bound using the current position */
1408   F = find_tension(state_list, env, x, 0, 1); /* BUG. repeated calculation */
1409   //printf("Start with x = %g (F = %g)\n", x, F);
1410   k = accel_k(transition->k, F, env, transition->k_params);
1411   //printf("x %g\tF %g\tk %g\n", x, F, k);
1412   dtU = P1(max_prob, N_INIT(transition)) / k;  /* upper bound on dt */
1413   if (dtU > max_dt) {
1414     //printf("overshot max_dt\n");
1415     dtU = max_dt;
1416   }
1417   if (v == 0) /* discrete stepping, so force is temporarily constant */
1418     return dtU;
1419
1420   /* set a lower bound on dt too */
1421   dtL = 0.0;
1422
1423   /* The dt determined above may produce illegitimate forces or ks.
1424    * Reduce the upper bound until we have valid ks. */
1425   dt = dtU;
1426   F = find_tension(state_list, env, x+v*dt, 0, 1);
1427   while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1428     dtU /= 2.0;
1429     dt = dtU;
1430     F = find_tension(state_list, env, x+v*dt, 0, 1);
1431   }
1432   //printf("Try for dt = %g (F = %g)\n", dt, F);
1433   k = accel_k(transition->k, F, env, transition->k_params);
1434   /* returned k may be -1.0 */
1435   //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1436   while (k == -1.0) { /* reduce step until we hit a valid k */
1437     dtU /= 2.0;
1438     dt = dtU; /* hopefully, we can use the max dt, see if it works */
1439     F = find_tension(state_list, env, x+v*dt, 0, 1);
1440     //printf("Try for dt = %g (F = %g)\n", dt, F);
1441     k = accel_k(transition->k, F, env, transition->k_params);
1442     //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1443   }
1444   assert(dtU > 1e-14);      /* timestep to valid k too small */
1445   /* safe timestep back from x+dtU */
1446   dtUCur = P1(max_prob, N_INIT(transition)) / k;
1447   if (dtUCur >= dt)
1448     return dt; /* dtU is safe. */
1449
1450   /* dtU wasn't safe, lets see what would be. */
1451   while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1452     dt = (dtU + dtL) / 2.0;
1453     F = find_tension(state_list, env, x+v*dt, 0, 1);
1454     //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1455     k = accel_k(transition->k, F, env, transition->k_params);
1456     dtCur = P1(max_prob, N_INIT(transition)) / k;
1457     //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1458     if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1459       dtL = dt;
1460     else if (dtCur < dt) {  /* unsafe timestep back from x+dt, but...    */
1461       dtU = dt;             /* ... stepping out only dtCur would be safe */
1462       dtUCur = dtCur;
1463     } else
1464       break; /* dtCur = dt */
1465   }
1466   return MAX(dtUCur, dtL);
1467 }
1468
1469 @
1470
1471 Determine $dt$ for an array of potentially different folded domains.
1472 We apply [[search_dt()]] to each populated domain to find the maximum
1473 $dt$ the domain can handle.  The final $dt$ is less than the
1474 individual domain $dt$s (to satisfy [[search_dt()]], see above).  If
1475 $v>0$ (continuous pulling), $dt$ is also reduced to satisfy
1476 $|F(x+v*dt)-F(x)|<dF_\text{max}$, which limits errors due to our
1477 transition rate assumption that the force does not change appreciably
1478 over a single timestep.
1479 <<determine dt>>=
1480 <<search dt>>
1481 double determine_dt(list_t *state_list, list_t *transition_list,
1482                     environment_t *env, double F, double x, double max_dF,
1483                     double max_prob, double dt_max, double v, int transition)
1484 { /* Returns the timestep dt in seconds.
1485    * Takes the state and transition lists, a pointer to the environment,
1486    * the force F(x) in N, an extension x in m, max_dF in N,
1487    * max_prob between 0 and 1, dt_max in s, v in m/s, and transition in {0,1}*/
1488   double dt=dt_max, new_dt, F_new;
1489   assert(max_dF > 0.0);
1490   assert(max_prob > 0.0); assert(max_prob < 1.0);
1491   assert(dt_max > 0.0);
1492   assert(state_list != NULL);
1493   assert(transition_list != NULL);
1494   transition_list = head(transition_list);
1495
1496   if (v > 0) {
1497     /* Ensure |dF| < max_dF */
1498     F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1499     while (fabs(F_new - F) > max_dF) {
1500       dt /= 2.0;
1501       F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1502     }
1503   }
1504
1505   /* Ensure all individual domains pass search_dt() */
1506   while (transition_list != NULL) {
1507     if (T(transition_list)->initial_state->num_domains > 0) {
1508       new_dt = search_dt(T(transition_list),state_list,env,x,max_prob,dt,v);
1509       dt = MIN(dt, new_dt);
1510     }
1511     transition_list = transition_list->next;
1512   }
1513   return dt;
1514 }
1515
1516 @
1517
1518 \subsection{State data}
1519 \label{sec.state_data}
1520
1521 The domains exist in one of an array of possible states.  Each state
1522 has a tension model which determines it's force/extension curve
1523 (possibly [[null]] for rigid states, see Appendix
1524 \ref{sec.null_tension}).
1525
1526 Groups are, for example, for two domain states with WLC tensions; one
1527 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1528 $L=28\U{nm}$.  Since the persistence lengths are the same, you would
1529 like to add the contour lengths of all the domains in \emph{both}
1530 states, and plug that total contour length into a single WLC
1531 evaluation (see Section \ref{sec.find_tension}).
1532 <<state definitions>>=
1533 typedef struct state_struct {
1534   char *name;                    /* name string                           */
1535   one_dim_fn_t *tension_handler; /* tension handler                       */
1536   one_dim_fn_t *inverse_tension_handler; /* inverse tension handler       */
1537   int tension_group;             /* for grouping similar tension states   */
1538   void *tension_params;          /* pointer to tension parameters         */
1539   destroy_data_func_t *destroy_tension;
1540   int num_domains;               /* number of domains currently in state  */
1541   /* cached lookups generated from the above data */
1542   int num_out_trans;             /* length of out_trans array             */
1543   struct transition_struct **out_trans; /* transitions leaving this state */
1544   int num_grouped_states;        /* length of grouped states array        */
1545   struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1546 } state_t;
1547
1548 /* get the state data for the current list node */
1549 #define S(list) ((state_t *)(list)->d)
1550
1551 @ [[tension_params]] is a pointer to the parameters used by the
1552 handler function when determining the tension.  [[destroy_tension]]
1553 points to a function for freeing [[tension_params]].  We it in the
1554 state structure so that [[destroy_state]] and will have an easier job
1555 cleaning up.
1556
1557 [[create_]] and [[destroy_state]] are simple wrappers around
1558 [[malloc]] and [[free]].
1559 <<create/destroy state>>=
1560 state_t *create_state(char *name,
1561                       one_dim_fn_t *tension_handler,
1562                       one_dim_fn_t *inverse_tension_handler,
1563                       int tension_group,
1564                       void *tension_params,
1565                       destroy_data_func_t *destroy_tension,
1566                       int num_domains)
1567 {
1568   state_t *ret = (state_t *)malloc(sizeof(state_t));
1569   assert(ret != NULL);
1570   /* make a copy of the name, so we aren't dependent on the original */
1571   ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1572   assert(ret->name != NULL);
1573   strcpy(ret->name, name); /* we know ret->name is long enough */
1574
1575   ret->tension_handler = tension_handler;
1576   ret->inverse_tension_handler = inverse_tension_handler;
1577   ret->tension_group = tension_group;
1578   ret->tension_params = tension_params;
1579   ret->destroy_tension = destroy_tension;
1580   ret->num_domains = num_domains;
1581
1582   ret->num_out_trans = 0;
1583   ret->out_trans = NULL;
1584   ret->num_grouped_states = 0;
1585   ret->grouped_states = NULL;
1586
1587 #ifdef DEBUG
1588   fprintf(stderr, "generated state %s (%p) with tension handler %p, inverse handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->inverse_tension_handler, ret->tension_params, ret->tension_group);
1589 #endif
1590   return ret;
1591 }
1592
1593 void destroy_state(state_t *state)
1594 {
1595   if (state) {
1596 #ifdef DEBUG
1597     fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1598 #endif
1599     if (state->name)
1600       free(state->name);
1601     if (state->destroy_tension)
1602       (*state->destroy_tension)(state->tension_params);
1603     if (state->out_trans)
1604       free(state->out_trans);
1605     if (state->grouped_states)
1606       free(state->grouped_states);
1607     free(state);
1608   }
1609 }
1610
1611 @
1612
1613 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1614 so we also define a convenience function for cleaning up.
1615 <<destroy state list>>=
1616 void destroy_state_list(list_t *state_list)
1617 {
1618   state_list = head(state_list);
1619   while (state_list != NULL)
1620     destroy_state((state_t *) pop(&state_list));
1621 }
1622
1623 @
1624
1625 We can also define a convenient way to get a state index from it's name.
1626 <<state index>>=
1627 int state_index(list_t *state_list, char *name)
1628 {
1629   int i=0;
1630   state_list = head(state_list);
1631   while (state_list != NULL) {
1632     if (STRMATCH(S(state_list)->name, name)) return i;
1633     state_list = state_list->next;
1634     i++;
1635   }
1636   return -1;
1637 }
1638
1639 @
1640
1641 \subsection{Transition data}
1642 \label{sec.transition_data}
1643
1644 Once you've created states (Sections \ref{sec.main},
1645 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1646 create transitions between them.
1647 <<transition definitions>>=
1648 typedef struct transition_struct {
1649   state_t *initial_state;
1650   state_t *final_state;
1651   k_func_t *k;            /* transition rate model   */
1652   void *k_params;         /* pointer to k parameters */
1653   destroy_data_func_t *destroy_k;
1654 } transition_t;
1655
1656 /* get the transition data for the current list node */
1657 #define T(list) ((transition_t *)(list)->d)
1658
1659 /* get the number of initial state population for a transition state */
1660 #define N_INIT(transition) ((transition)->initial_state->num_domains)
1661
1662 <<PN definition>>
1663 <<P1 definition>>
1664 @ [[k]] is a pointer to the function determining the transition rate
1665 for a given tension.  [[k_params]] and [[destroy_k]] have similar
1666 roles with regards to [[k]] as [[tension_params]] and
1667 [[destroy_tension]] have with regards to [[tension_handler]] (see
1668 Section \ref{sec.state_data}).
1669
1670 [[create_]] and [[destroy_transition]] are simple wrappers around
1671 [[malloc]] and [[free]].
1672 <<create/destroy transition>>=
1673 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1674                                 k_func_t *k,
1675                                 void *k_params,
1676                                 destroy_data_func_t *destroy_k)
1677 {
1678   transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1679   assert(ret != NULL);
1680   assert(initial_state != NULL);
1681   assert(final_state != NULL);
1682   ret->initial_state = initial_state;
1683   ret->final_state = final_state;
1684   ret->k = k;
1685   ret->k_params = k_params;
1686   ret->destroy_k = destroy_k;
1687 #ifdef DEBUG
1688   fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1689 #endif
1690   return ret;
1691 }
1692
1693 void destroy_transition(transition_t *transition)
1694 {
1695   if (transition) {
1696 #ifdef DEBUG
1697     fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1698 #endif
1699     if (transition->destroy_k)
1700       (*transition->destroy_k)(transition->k_params);
1701     free(transition);
1702   }
1703 }
1704
1705 @
1706
1707 We'll be storing the transitions in a list (see Appendix
1708 \ref{sec.lists}), so we also define a convenience function for
1709 cleaning up.
1710 <<destroy transition list>>=
1711 void destroy_transition_list(list_t *transition_list)
1712 {
1713   transition_list = head(transition_list);
1714   while (transition_list != NULL)
1715     destroy_transition((transition_t *) pop(&transition_list));
1716 }
1717
1718 @
1719
1720 \subsection{Graphviz output}
1721 \label{sec.graphviz_output}
1722
1723 <<print state graph>>=
1724 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1725 {
1726   state_list = head(state_list);
1727   transition_list = head(transition_list);
1728   fprintf(file, "digraph states {\n  node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1729   while (1) {
1730     fprintf(file, "  n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1731     if (state_list->next == NULL) break;
1732     state_list = state_list->next;
1733   }
1734   fprintf(file, "\n");
1735   while (transition_list != NULL) {
1736     fprintf(file, "  n%d -> n%d;\n", state_index(state_list, T(transition_list)->initial_state->name), state_index(state_list, T(transition_list)->final_state->name));
1737     transition_list = transition_list->next;
1738   }
1739   fprintf(file, "}\n");
1740 }
1741 @
1742
1743 \subsection{Domain model and group handling}
1744
1745 <<group functions>>=
1746 <<crosslink groups>>
1747 <<get active groups>>
1748 @
1749
1750 Scan through a list of states and construct an array of all the
1751 active groups.
1752 <<get active groups>>=
1753 void get_active_groups(list_t *state_list,
1754                        int *pNum_active_groups,
1755                        one_dim_fn_t ***pPer_group_handlers,
1756                        one_dim_fn_t ***pPer_group_inverse_handlers,
1757                        void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1758 {
1759   void **active_groups=NULL;
1760   one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1761   one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1762   void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1763   int i, j, num_domains, group, old_group, num_active_groups=0;
1764   list_t *active_states_list=NULL;
1765   tension_handler_data_t *tdata=NULL;
1766   state_t *state;
1767
1768   /* get all the active groups in a list */
1769   state_list = head(state_list);
1770 #ifdef DEBUG
1771   fprintf(stderr, "%s:\t", __FUNCTION__);
1772   list_print(stderr, state_list, "states");
1773 #endif
1774   while (state_list != NULL) {
1775     state = S(state_list);
1776     handler = state->tension_handler;
1777     inverse_handler = state->inverse_tension_handler;
1778     group = state->tension_group;
1779     num_domains = state->num_domains;
1780     if (list_index(active_states_list, state) == -1) {
1781       /* we haven't added this state (or it's associated group) yet */
1782       if (num_domains > 0 && handler != NULL) { /* add the state */
1783         num_active_groups++;
1784         push(&active_states_list, state, 1);
1785 #ifdef DEBUG
1786         fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1787 #endif
1788         for (i=0; i < state->num_grouped_states; i++) {
1789           if (state->grouped_states[i]->num_domains > 0) {
1790             /* add this related (and active) state now */
1791             assert(state->grouped_states[i]->tension_handler == handler);
1792             assert(state->grouped_states[i]->tension_group == group);
1793             push(&active_states_list, state->grouped_states[i], 1);
1794 #ifdef DEBUG
1795             fprintf(stderr, "%s:\tactive grouped %s state (%p) with %d domain(s)\n", __FUNCTION__, state->grouped_states[i]->name, state->grouped_states[i], state->grouped_states[i]->num_domains);
1796 #endif
1797           }
1798         }
1799       } /* else empty state or NULL handler */
1800     } /* else state already added */
1801     state_list = state_list->next;
1802   }
1803 #ifdef DEBUG
1804   fprintf(stderr, "%s:\t", __FUNCTION__);
1805   list_print(stderr, active_states_list, "active states");
1806 #endif
1807
1808   assert(num_active_groups <= list_length(active_states_list));
1809   assert(num_active_groups > 0);
1810
1811   /* allocate the arrays we'll be returning */
1812   per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1813   assert(per_group_handlers != NULL);
1814   per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1815   assert(per_group_inverse_handlers != NULL);
1816   per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1817   assert(per_group_data != NULL);
1818 #ifdef DEBUG
1819   fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1820 #endif
1821   for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1822     per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1823     assert(per_group_data[i] != NULL);
1824 #ifdef DEBUG
1825     fprintf(stderr, " %d,%p", i, per_group_data[i]);
1826 #endif
1827   }
1828 #ifdef DEBUG
1829     fprintf(stderr, "\n");
1830 #endif
1831
1832   /* populate the arrays we'll be returning */
1833   active_states_list = head(active_states_list);
1834   old_handler = NULL;
1835   old_inverse_handler = NULL;
1836   j = -1; /* j is the current active _group_ index */
1837   while (active_states_list != NULL) {
1838     state = (state_t *) pop(&active_states_list);
1839     handler = state->tension_handler;
1840     inverse_handler = state->inverse_tension_handler;
1841     group = state->tension_group;
1842     num_domains = state->num_domains;
1843     if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1844       j++; /* move to the next group */
1845       tdata = (tension_handler_data_t *) per_group_data[j];
1846       per_group_handlers[j] = handler;
1847       per_group_inverse_handlers[j] = inverse_handler;
1848       tdata->group_tension_params = NULL;
1849       tdata->env = NULL;
1850       tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1851     }
1852 #ifdef DEBUG
1853     fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, inverse_handler, group, num_domains);
1854 #endif
1855     for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1856       push(&tdata->group_tension_params, state->tension_params, 1);
1857     }
1858 #ifdef DEBUG
1859     fprintf(stderr, "%s:\tpushed %d copies of %p onto data %p's param list %p\n", __FUNCTION__, num_domains, state->tension_params, tdata, tdata->group_tension_params);
1860 #endif
1861     old_handler = handler;
1862     old_inverse_handler = inverse_handler;
1863     old_group = group;
1864   }
1865
1866   /* free old groups */
1867   if (*pPer_group_handlers != NULL)
1868     free(*pPer_group_handlers);
1869   if (*pPer_group_inverse_handlers != NULL)
1870     free(*pPer_group_inverse_handlers);
1871   if (*pPer_group_data != NULL) {
1872     for (i=0; i<*pNum_active_groups; i++)
1873       free((*pPer_group_data)[i]);
1874     free(*pPer_group_data);
1875   }
1876
1877   /* set new groups */
1878 #ifdef DEBUG
1879   fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
1880 #endif
1881   *pNum_active_groups = num_active_groups;
1882   *pPer_group_handlers = per_group_handlers;
1883   *pPer_group_inverse_handlers = per_group_inverse_handlers;
1884   *pPer_group_data = per_group_data;
1885 }
1886 @
1887
1888 <<crosslink groups>>=
1889 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1890 {
1891   list_t *list, *out_trans=NULL, *associates=NULL;
1892   one_dim_fn_t *handler;
1893   int i, n, group;
1894
1895   state_groups = head(state_groups);
1896   while (state_groups != NULL) {
1897     /* find transitions leaving this state */
1898 #ifdef DEBUG
1899     fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1900 #endif
1901     list = head(transition_list);
1902     while (list != NULL) {
1903       if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1904         push(&out_trans, T(list), 1);
1905       }
1906       list = list->next;
1907     }
1908     n = list_length(out_trans);
1909     S(state_groups)->num_out_trans = n;
1910     S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1911     assert(S(state_groups)->out_trans != NULL);
1912     for (i=0; i<n; i++) {
1913       S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1914     }
1915
1916     /* find states grouped with this state */
1917 #ifdef DEBUG
1918     fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1919 #endif
1920     handler = S(state_groups)->tension_handler;
1921     group = S(state_groups)->tension_group;
1922     list = head(state_groups);
1923     while (list != NULL) {
1924       if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1925         push(&associates, S(list), 1);
1926       }
1927       list = list->next;
1928     }
1929     n = list_length(associates);
1930     S(state_groups)->num_grouped_states = n;
1931     S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1932     assert(S(state_groups)->grouped_states != NULL);
1933     for (i=0; i<n; i++) {
1934       S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1935     }
1936   state_groups = state_groups->next;
1937   }
1938 }
1939
1940 @
1941
1942 \section{String parsing}
1943
1944 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1945 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1946 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1947 need to take care of parsing those parameters themselves.
1948 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1949
1950 <<parse.h>>=
1951 #ifndef PARSE
1952 #define PARSE
1953 <<license comment>>
1954 <<parse definitions>>
1955 <<parse declarations>>
1956 #endif /* PARSE */
1957 @
1958
1959 <<parse module makefile lines>>=
1960 NW_SAWSIM_MODS += parse
1961 CHECK_BINS += check_parse
1962 check_parse_MODS =
1963 @
1964
1965 <<parse definitions>>=
1966 #define SEP ',' /* argument separator character */
1967 @
1968
1969 <<parse declarations>>=
1970 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1971                        int *num, char ***string_array);
1972 extern void free_parsed_list(int num, char **string_array);
1973 @
1974
1975 [[parse_list_string]] allocates memory, don't forget to free it
1976 afterward with [[free_parsed_list]].  It does not alter the original.
1977
1978 The string may start off with a [[deeper]] character
1979 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1980 [[x,y]], where the pointer is one character in on the copied string.
1981 However, when we go to free the memory, we need a pointer to the
1982 beginning of the string.  In order to accommodate this for a string
1983 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1984 the first $N$ elements point to the separated arguments, and let the
1985 last element point to the start of the copied string regardless of
1986 braces.
1987 <<parse delimited list string functions>>=
1988 /* TODO, split out into parse.hc */
1989 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1990 {
1991   int i=0, depth = 0;
1992   while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1993     if (string[i] == deeper) {depth++;}
1994     else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1995     i++;
1996   }
1997   return i;
1998 }
1999
2000 void parse_list_string(char *string, char sep, char deeper, char shallower,
2001                        int *num, char ***string_array)
2002 {
2003   char *str=NULL, **ret=NULL;
2004   int i, j, n;
2005   if (string==NULL || strlen(string) == 0) {    /* handle the trivial cases */
2006     *num = 0;
2007     *string_array = NULL;
2008     return;
2009   }
2010   /* make a copy of the string, so we don't change the original */
2011   str = (char *)malloc(sizeof(char)*(strlen(string)+1));
2012   assert(str != NULL);
2013   strcpy(str, string); /* we know str is long enough */
2014   /* count the number of regions, so we can allocate pointers to them */
2015   i=-1; n=0;
2016   do {
2017     n++; i++; /* move on to next argument */
2018     i += next_delim_index(str+i, sep, deeper, shallower);
2019     //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
2020   } while (str[i] != '\0');
2021   ret = (char **)malloc(sizeof(char *)*(n+1));
2022   assert(ret != NULL);
2023   /* replace the separators with '\0' & assign pointers */
2024   ret[n] = str; /* point to the front of the copied string */
2025   j=0;
2026   ret[0] = str;
2027   for(i=1; i<n; i++) {
2028     j += next_delim_index(str+j, sep, deeper, shallower);
2029     str[j++] = '\0';
2030     ret[i] = str+j; /* point to the separated arguments */
2031   }
2032   /* strip off bounding braces, if any */
2033   for(i=0; i<n; i++) {
2034     if (ret[i][0] == deeper) {
2035       ret[i][0] = '\0';
2036       ret[i]++;
2037     }
2038     if (ret[i][strlen(ret[i])-1] == shallower)
2039       ret[i][strlen(ret[i])-1] = '\0';
2040   }
2041   *num = n;
2042   *string_array = ret;
2043 }
2044
2045 void free_parsed_list(int num, char **string_array)
2046 {
2047   if (string_array) {
2048     /* frees all items, since they were allocated as a single string*/
2049     free(string_array[num]);
2050     /* and free the array of pointers */
2051     free(string_array);
2052   }
2053 }
2054 @
2055
2056 \subsection{Parsing implementation}
2057
2058 <<parse.c>>=
2059 <<license comment>>
2060 <<parse includes>>
2061 #include "parse.h"
2062 <<parse delimited list string functions>>
2063 @
2064
2065 <<parse includes>>=
2066 #include <assert.h> /* assert()                */
2067 #include <stdlib.h> /* NULL                    */
2068 #include <stdio.h>  /* fprintf(), stdout       *//*!!*/
2069 #include <string.h> /* strlen()                */
2070 #include "parse.h"
2071 @
2072
2073 \subsection{Parsing unit tests}
2074
2075 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2076 <<check-parse.c>>=
2077 <<license comment>>
2078 <<parse check includes>>
2079 <<string comparison definition and globals>>
2080 <<check relative error>>
2081 <<parse test suite>>
2082 <<main check program>>
2083 @
2084
2085 <<parse check includes>>=
2086 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2087 #include <stdio.h>  /* printf()                      */
2088 #include <assert.h> /* assert()                      */
2089 #include <string.h> /* strlen()                      */
2090 <<check includes>>
2091 #include "parse.h"
2092 @
2093
2094 <<parse test suite>>=
2095 <<parse list string tests>>
2096 <<parse suite definition>>
2097 @
2098
2099 <<parse suite definition>>=
2100 Suite *test_suite (void)
2101 {
2102   Suite *s = suite_create ("parse");
2103   <<parse list string test case defs>>
2104
2105   <<parse list string test case adds>>
2106   return s;
2107 }
2108 @
2109
2110 <<parse list string tests>>=
2111
2112 /*
2113 START_TEST(test_next_delim_index)
2114 {
2115   fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2116   fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2117   fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2118   fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2119   fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2120 }
2121 END_TEST
2122 */
2123 START_TEST(test_parse_list_null)
2124 {
2125   int num_param_args;
2126   char **param_args;
2127   parse_list_string(NULL, SEP, '{', '}',
2128                     &num_param_args, &param_args);
2129   fail_unless(num_param_args == 0, NULL);
2130   fail_unless(param_args == NULL, NULL);
2131 }
2132 END_TEST
2133 START_TEST(test_parse_list_single_simple)
2134 {
2135   int num_param_args;
2136   char **param_args;
2137   parse_list_string("arg", SEP, '{', '}',
2138                     &num_param_args, &param_args);
2139   fail_unless(num_param_args == 1, NULL);
2140   fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2141 }
2142 END_TEST
2143 START_TEST(test_parse_list_single_compound)
2144 {
2145   int num_param_args;
2146   char **param_args;
2147   parse_list_string("{x,y,z}", SEP, '{', '}',
2148                     &num_param_args, &param_args);
2149   fail_unless(num_param_args == 1, NULL);
2150   fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2151 }
2152 END_TEST
2153 START_TEST(test_parse_list_double_simple)
2154 {
2155   int num_param_args;
2156   char **param_args;
2157   parse_list_string("abc,def", SEP, '{', '}',
2158                     &num_param_args, &param_args);
2159   fail_unless(num_param_args == 2, NULL);
2160   fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2161   fail_unless(STRMATCH(param_args[1],"def"), NULL);
2162 }
2163 END_TEST
2164 START_TEST(test_parse_list_compound)
2165 {
2166   int num_param_args;
2167   char **param_args;
2168   parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2169                     &num_param_args, &param_args);
2170   fail_unless(num_param_args == 2, NULL);
2171   fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2172   fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2173 }
2174 END_TEST
2175 @
2176 <<parse list string test case defs>>=
2177 TCase *tc_parse_list_string = tcase_create("parse list string");
2178 @
2179 <<parse list string test case adds>>=
2180 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2181 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2182 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2183 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2184 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2185 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2186 suite_add_tcase(s, tc_parse_list_string);
2187 @
2188
2189
2190 \section{Unit tests}
2191
2192 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2193 <<checks>>=
2194 <<includes>>
2195 <<check includes>>
2196 <<definitions>>
2197 <<globals>>
2198 <<check globals>>
2199 <<check relative error>>
2200 <<functions>>
2201 <<test suite>>
2202 <<main check program>>
2203 @
2204
2205 <<check includes>>=
2206 #include <check.h>
2207 @
2208
2209 <<check globals>>=
2210 @
2211
2212 <<test suite>>=
2213 <<determine dt tests>>
2214 <<happens tests>>
2215 <<does domain transition tests>>
2216 <<randomly unfold domains tests>>
2217 <<suite definition>>
2218 @
2219
2220 <<suite definition>>=
2221 Suite *test_suite (void)
2222 {
2223   Suite *s = suite_create ("sawsim");
2224   <<determine dt test case defs>>
2225   <<happens test case defs>>
2226   <<does domain transition test case defs>>
2227   <<randomly unfold domains test case defs>>
2228
2229   <<determine dt test case adds>>
2230   <<happens test case adds>>
2231   <<does domain transition test case adds>>
2232   <<randomly unfold domains test case adds>>
2233
2234 /*
2235   tcase_add_checked_fixture(tc_strip_address,
2236                             setup_strip_address,
2237                             teardown_strip_address);
2238 */
2239   return s;
2240 }
2241 @
2242
2243 <<main check program>>=
2244 int main(void)
2245 {
2246   int number_failed;
2247   Suite *s = test_suite();
2248   SRunner *sr = srunner_create(s);
2249   srunner_run_all(sr, CK_ENV);
2250   number_failed = srunner_ntests_failed(sr);
2251   srunner_free(sr);
2252   return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2253 }
2254 @
2255
2256
2257 \subsection{Model tests}
2258
2259 Check the searching with [[linear_k]].
2260 Check overwhelming force treatment with the heavyside-esque step [[?]].
2261
2262 Hmm, I'm not really sure what I was doing here...
2263 <<determine dt tests>>=
2264 double linear_k(double F, environment_t *env, void *params)
2265 {
2266   double Fnot = *(double *)params;
2267   return Fnot+F;
2268 }
2269
2270 START_TEST(test_determine_dt_linear_k)
2271 {
2272   state_t folded={0};
2273   transition_t unfold={0};
2274   environment_t env = {0};
2275   double F[]={0,1,2,3,4,5,6};
2276   double dt_max=3.0, Fnot=3.0, x=1.0, max_p=0.1, max_dF=0.1, v=1.0;
2277   list_t *state_list=NULL, *transition_list=NULL;
2278   int i;
2279   env.T = 300.0;
2280   folded.tension_handler = &hooke_handler;
2281   folded.num_domains = 1;
2282   unfold.initial_state = &folded;
2283   unfold.k = linear_k;
2284   unfold.k_params = &Fnot;
2285   push(&state_list, &folded, 1);
2286   push(&transition_list, &unfold, 1);
2287   for( i=0; i < sizeof(F)/sizeof(double); i++) {
2288     //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL);
2289   }
2290 }
2291 END_TEST
2292 @
2293 <<determine dt test case defs>>=
2294 TCase *tc_determine_dt = tcase_create("Determine dt");
2295 @
2296 <<determine dt test case adds>>=
2297 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2298 suite_add_tcase(s, tc_determine_dt);
2299 @
2300
2301 <<happens tests>>=
2302 @
2303 <<happens test case defs>>=
2304 @
2305 <<happens test case adds>>=
2306 @
2307
2308 <<does domain transition tests>>=
2309 @
2310 <<does domain transition test case defs>>=
2311 @
2312 <<does domain transition test case adds>>=
2313 @
2314
2315 <<randomly unfold domains tests>>=
2316 @
2317 <<randomly unfold domains test case defs>>=
2318 @
2319 <<randomly unfold domains test case adds>>=
2320 @
2321
2322
2323 \section{Balancing group extensions}
2324 \label{sec.tension_balance}
2325
2326 For a given total extension $x$ (of the piezo), the various domain
2327 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2328 amounts, and we need to tweak the portion that each extends to balance
2329 the tension amongst the active groups.  Since the tension balancing is
2330 separable from the bulk of the simulation, we break it out into a
2331 separate module.  The interface is defined in [[tension_balance.h]],
2332 the implementation in [[tension_balance.c]], and the unit testing in
2333 [[check_tension_balance.c]]
2334
2335 <<tension-balance.h>>=
2336 #ifndef TENSION_BALANCE_H
2337 #define TENSION_BALANCE_H
2338 <<license comment>>
2339 <<tension balance function declaration>>
2340 #endif /* TENSION_BALANCE_H */
2341 @
2342
2343 <<tension balance functions>>=
2344 <<one dimensional bracket>>
2345 <<one dimensional solve>>
2346 <<x of x0>>
2347 <<group stiffness function>>
2348 <<chain stiffness function>>
2349 <<full chain stiffness function>>
2350 <<tension balance function>>
2351 @
2352
2353 <<tension balance module makefile lines>>=
2354 NW_SAWSIM_MODS += tension_balance
2355 CHECK_BINS += check_tension_balance
2356 check_tension_balance_MODS =
2357 @
2358
2359 The entire force balancing problem reduces to a solving two nested
2360 one-dimensional root-finding problems.  First we define the one
2361 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2362 populated group).  $x(x_0)$ is also strictly monotonically increasing,
2363 so we can solve for $x_0$ such that $\sum_i x_i = x$.  [[last_x]]
2364 stores the last successful value of $x$, since we expect to be taking
2365 small steps ($x \approx$ [[last_x]]).  Setting [[last_x = -1]]
2366 indicates that the groups have changed.
2367 <<tension balance function declaration>>=
2368 double tension_balance(int num_tension_groups,
2369                        one_dim_fn_t **tension_handlers,
2370                        one_dim_fn_t **inverse_tension_handlers,
2371                        void **params, /* array of pointers to tension_handler_data_t */
2372                        double last_x,
2373                        double x,
2374                        double *stiffness);
2375 @
2376 <<tension balance function>>=
2377 double tension_balance(int num_tension_groups,
2378                        one_dim_fn_t **tension_handlers,
2379                        one_dim_fn_t **inverse_tension_handlers,
2380                        void **params, /* array of pointers to tension_handler_data_t */
2381                        double last_x,
2382                        double x,
2383                        double *stiffness)
2384 {
2385   static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2386   double F, xo_guess, xo, lb, ub;
2387   double min_relx=1e-6, min_rely=1e-6;
2388   int max_steps=200, i;
2389
2390   assert(num_tension_groups > 0);
2391   assert(tension_handlers != NULL);
2392   assert(params != NULL);
2393   assert(num_tension_groups > 0);
2394
2395   if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2396     /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2397     if (x_xo_data.xi != NULL)
2398       free(x_xo_data.xi);
2399     /* construct new x_xo_data */
2400     x_xo_data.n_groups = num_tension_groups;
2401     x_xo_data.tension_handler = tension_handlers;
2402     x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2403     x_xo_data.handler_data = params;
2404 #ifdef DEBUG
2405     fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p), x_of_xo_data %p\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0], &x_xo_data);
2406 #endif
2407     x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2408     assert(x_xo_data.xi != NULL);
2409     for (i=0; i<num_tension_groups; i++)
2410       x_xo_data.xi[i] = last_x;
2411   }
2412
2413   if (num_tension_groups == 1) { /* only one group, no balancing required */
2414     xo = x;
2415     x_xo_data.xi[0] = xo;
2416 #ifdef DEBUG
2417     fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2418 #endif
2419   } else {
2420 #ifdef DEBUG
2421     fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2422     for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2423     fprintf(stderr, "\n");
2424 #endif
2425     if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2426       if (x == 0) xo_guess = length_scale;
2427       else        xo_guess = x/num_tension_groups;
2428 #ifdef DEBUG
2429       fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2430 #endif
2431       oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2432     } else { /* work off the last balanced point */
2433 #ifdef DEBUG
2434       fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2435 #endif
2436       if (x > last_x) {
2437         lb = x_xo_data.xi[0];
2438         ub = x_xo_data.xi[0]+ x-last_x;   /* apply all change to x0 */
2439       } else if (x < last_x) {
2440         lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2441         ub = x_xo_data.xi[0];
2442       } else { /* x == last_x */
2443 #ifdef DEBUG
2444         fprintf(stderr, "not moving\n");
2445 #endif
2446         lb = x_xo_data.xi[0];
2447         ub = x_xo_data.xi[0];
2448       }
2449     }
2450 #ifdef DEBUG
2451     fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2452             __FUNCTION__, x, lb, ub);
2453 #endif
2454     xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2455     F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2456 #ifdef DEBUG
2457     if (num_tension_groups > 1) {
2458       fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2459       for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2460       fprintf(stderr, "\n");
2461     }
2462 #endif
2463   }
2464
2465   F = (*tension_handlers[0])(xo, params[0]);
2466
2467   if (stiffness != NULL) {
2468     *stiffness = chain_stiffness(&x_xo_data, min_relx);
2469     if (*stiffness == 0.0) { /* re-calculate based on full chain */
2470       *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
2471                                         inverse_tension_handlers, params,
2472                                         last_x, x, min_relx, F);
2473     }
2474   }
2475   return F;
2476 }
2477
2478 @
2479
2480 <<tension balance internal definitions>>=
2481 <<x of x0 definitions>>
2482 @
2483
2484 <<x of x0 definitions>>=
2485 typedef struct x_of_xo_data_struct {
2486   int n_groups;
2487   one_dim_fn_t **tension_handler;        /* array of fn pointers      */
2488   one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers      */
2489   void **handler_data; /* array of pointers to tension_handler_data_t */
2490   double *xi;                            /* array of group extensions */
2491 } x_of_xo_data_t;
2492 @
2493
2494 <<x of x0>>=
2495 double x_of_xo(double xo, void *pdata)
2496 {
2497   x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2498   double F, x=0, xi, xi_guess, lb, ub, slope;
2499   int i;
2500   double min_relx=1e-6, min_rely=1e-6;
2501   int max_steps=200;
2502   assert(data->n_groups > 0);
2503   data->xi[0] = xo;
2504 #ifdef DEBUG
2505   fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p (x_xo_data %p)\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0], data);
2506   fflush(stderr);
2507 #endif
2508   F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2509 #ifdef DEBUG
2510   fprintf(stderr, "%s:\tfinding x(x0=%g).  F0(x0) = %g\n", __FUNCTION__, xo, F);
2511   fflush(stderr);
2512 #endif
2513   x += xo;
2514   if (data->xi)  data->xi[0] = xo;
2515   for (i=1; i < data->n_groups; i++) {
2516     if (data->inverse_tension_handler[i] != NULL) {
2517       xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2518     } else { /* invert numerically */
2519       if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2520       else                   xi_guess = data->xi[i];
2521 #ifdef DEBUG
2522       fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]);
2523 #endif
2524       oneD_bracket(data->tension_handler[i], data->handler_data[i], F,  xi_guess, &lb, &ub);
2525       xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2526     }
2527     x += xi;
2528     data->xi[i] = xi;
2529 #ifdef DEBUG
2530     fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2531 #endif
2532   }
2533 #ifdef DEBUG
2534   fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2535 #endif
2536   return x;
2537 }
2538 @
2539
2540 Determine the stiffness of a single tension group by taking a small
2541 step [[dx]] back from the position [[x]] for which we want the
2542 stiffness.  The touchy part is determining a good [[dx]] and ensuring
2543 the whole interval is on [[x>=0]].
2544 <<group stiffness function>>=
2545 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2546 {
2547   double dx, xl, F, Fl, stiffness;
2548   assert(x >= 0);
2549   assert(relx > 0 && relx < 1);
2550   if (x == 0 || relx == 0) {
2551     dx = 1e-12; /* HACK, how to get length scale? */
2552     xl = x-dx;
2553     if (xl < 0) {
2554       xl = 0;
2555       x = xl+dx;
2556     }
2557   } else {
2558     dx = x*relx;
2559     xl = x-dx;
2560   }
2561   F = (*tension_handler)(x, handler_data);
2562   Fl = (*tension_handler)(xl, handler_data);
2563   stiffness = (F-Fl)/dx;
2564   assert(stiffness >= 0);
2565   return stiffness;
2566 }
2567 @
2568
2569 Attempt to calculate the chain stiffness by adding group stiffnesses
2570 as springs in series.  In the case where a group has undefined
2571 stiffness (e.g. piston tension, Section \ref{sec.piston_tension}),
2572 this algorithm breaks down.  In that case, [[tension_balance()]]
2573 (\ref{sec.tension_balance}) should use [[full_chain_stiffness()]]
2574 which uses the full chain's $F(x)$ rather than that of the individual
2575 domains'.  We attempt the series approach first because, lacking the
2576 numerical inversion inside [[tension_balance()]], it is both faster
2577 and more accurate than the full-chain derivative.
2578 <<chain stiffness function>>=
2579 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2580 {
2581   double stiffness, gstiff;
2582   int i;
2583   /* add group stiffnesses in series */
2584   for (i=0; i < x_xo_data->n_groups; i++) {
2585 #ifdef DEBUG
2586     fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx);
2587     fflush(stderr);
2588 #endif
2589     gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2590     if (gstiff == 0.0);
2591       return 0.0;
2592     stiffness += 1.0/gstiff;
2593   }
2594   stiffness = 1.0 / stiffness;
2595   return stiffness;
2596 }
2597
2598 @
2599
2600 Determine the chain stiffness using only [[tension_balance()]].  This
2601 works around difficulties with tension models that have undefined
2602 stiffness (see discussion for [[chain_stiffness()]] above).
2603 <<full chain stiffness function>>=
2604 double full_chain_stiffness(int num_tension_groups,
2605                        one_dim_fn_t **tension_handlers,
2606                        one_dim_fn_t **inverse_tension_handlers,
2607                        void **params, /* array of pointers to tension_handler_data_t */
2608                        double last_x,
2609                        double x,
2610                        double relx,
2611                        double F /* already evaluated F(x) */)
2612 {
2613   double dx, xl, Fl, stiffness;
2614
2615   assert(x >= 0);
2616   assert(relx > 0 && relx < 1);
2617
2618   /* Other option for dx that we ignore because of our poor tension_balance()
2619    * resolution (i.e. bad slopes if you zoom in too close):
2620    *     if (last_x != -1.0 && last_x != x)  // last_x limited
2621    *       dx fabs(x-last_x);
2622    *     else
2623    *       dx = HUGE_VAL;
2624    *     if (1==1) {                 // constant max-value limited
2625    *       dx_p = 1e-12;
2626    *       dx = MIN(dx, dx_p);
2627    *     }
2628    *     if (x != 0 && relx != 0) {  // relx limited
2629    *       dx_p = x*relx;
2630    *       dx = MIN(dx, dx_p);
2631    *     }
2632    * TODO, determine which of (min_relx, min_rely, max_steps) is actually
2633    * limiting tension_balance.
2634    */
2635   dx = 1e-12; /* HACK, how to get length scale? */
2636
2637   xl = x-dx;
2638   if (xl >= 0) {
2639     Fl = tension_balance(num_tension_groups, tension_handlers,
2640                          inverse_tension_handlers, params,
2641                          last_x, xl, NULL);
2642   } else {
2643     xl = x;
2644     Fl = F;
2645     x += dx;
2646     F = tension_balance(num_tension_groups, tension_handlers,
2647                        inverse_tension_handlers, params,
2648                        last_x, x, NULL);
2649   }
2650   stiffness = (F-Fl)/dx;
2651   assert(stiffness >= 0);
2652   return stiffness;
2653 }
2654
2655 @
2656
2657 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2658 number of steps for monotonically increasing $f(x)$.  Simple
2659 bisection, so it's robust and converges fairly quickly.  You can set
2660 the maximum number of steps to take, but if you have enough processor
2661 power, it's better to limit your precision with the relative limits
2662 [[min_relx]] and [[y]].  These ensure that the width of the bracket is
2663 small on the length scale of it's larger side.  Note that \emph{both}
2664 relative limits must be satisfied for the search to stop.
2665 <<one dimensional function definition>>=
2666 /* equivalent to gsl_func_t */
2667 typedef double one_dim_fn_t(double x, void *params);
2668 @
2669 <<one dimensional solve declaration>>=
2670 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2671                   double min_relx, double min_rely, int max_steps,
2672                   double *pYx);
2673 @
2674
2675 <<one dimensional solve>>=
2676 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2677 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2678                   double min_relx, double min_rely, int max_steps,
2679                   double *pYx)
2680 {
2681   double yx, ylb, yub, x;
2682   int i=0;
2683   assert(ub >= lb);
2684   ylb = (*f)(lb, params);
2685   yub = (*f)(ub, params);
2686
2687   /* check some simple cases */
2688   if (ylb == yub) {
2689     if (ylb != y)  return HUGE_VAL; /* error! f(x)=y not bracketed */
2690     else           return lb; /* any x on the range [lb, ub] would work */
2691   }
2692   if (ylb == y) { x=lb; yx=ylb; goto end; }
2693   if (yub == y) { x=ub; yx=yub; goto end; }
2694
2695 #ifdef DEBUG
2696   fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2697 #endif
2698   assert(ylb < y); assert(yub > y);
2699   for (i=0; i<max_steps; i++) {
2700     x = (lb+ub)/2.0;
2701     yx = (*f)(x, params);
2702 #ifdef DEBUG
2703   fprintf(stderr, "%s:\tbracket: lb %g, x %g, ub %g\tylb %g, yx %g, yub %g\t(y %g)\n", __FUNCTION__, lb, x, ub, ylb, yx, yub, y);
2704 #endif
2705     if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2706     else if (yx > y)  { ub=x; yub=yx; }
2707     else /* yx < y */ { lb=x; ylb=yx; }
2708   }
2709  end:
2710   if (pYx) *pYx = yx;
2711 #ifdef DEBUG
2712   fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2713 #endif
2714   return x;
2715 }
2716 @
2717
2718 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2719 Generate bracketing $x$ values through bisection or exponential growth.
2720 <<one dimensional bracket declaration>>=
2721 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2722 @
2723
2724 <<one dimensional bracket>>=
2725 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2726 {
2727   double yx, step, x=xguess;
2728   int i=0;
2729 #ifdef DEBUG
2730   fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2731 #endif
2732   yx = (*f)(x, params);
2733 #ifdef DEBUG
2734   fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2735 #endif
2736   if (yx > y) {
2737     assert(x > 0.0);
2738     *ub = x;
2739     *lb = 0;
2740 #ifdef DEBUG
2741     fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2742 #endif
2743   } else {
2744     *lb = x;
2745     if (x == 0) x = length_scale/2.0; /* HACK */
2746     while (yx < y) {
2747       x *= 2.0;
2748       yx = (*f)(x, params);
2749 #ifdef DEBUG
2750       fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2751 #endif
2752     }
2753     *ub = x;
2754 #ifdef DEBUG
2755     fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2756 #endif
2757   }
2758   //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2759 }
2760 @
2761
2762 \subsection{Balancing implementation}
2763
2764 <<tension-balance.c>>=
2765 //#define DEBUG
2766 <<license comment>>
2767 <<tension balance includes>>
2768 #include "tension_balance.h"
2769
2770 double length_scale = 1e-10; /* HACK */
2771
2772 <<tension balance internal definitions>>
2773 <<tension balance functions>>
2774 @
2775
2776 <<tension balance includes>>=
2777 #include <assert.h> /* assert()          */
2778 #include <stdlib.h> /* NULL              */
2779 #include <math.h>   /* HUGE_VAL, macro constant, so don't need to link to math */
2780 #include <stdio.h>  /* fprintf(), stdout */
2781 #include "global.h"
2782 @
2783
2784 \subsection{Balancing unit tests}
2785
2786 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2787 <<check-tension-balance.c>>=
2788 <<license comment>>
2789 <<tension balance check includes>>
2790 <<tension balance test suite>>
2791 <<main check program>>
2792 @
2793
2794 <<tension balance check includes>>=
2795 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2796 #include <assert.h> /* assert()                      */
2797 <<check includes>>
2798 #include "global.h"
2799 #include "tension_balance.h"
2800 @
2801
2802 <<tension balance test suite>>=
2803 <<tension balance function tests>>
2804 <<tension balance suite definition>>
2805 @
2806
2807 <<tension balance suite definition>>=
2808 Suite *test_suite (void)
2809 {
2810   Suite *s = suite_create ("tension balance");
2811   <<tension balance function test case defs>>
2812
2813   <<tension balance function test case adds>>
2814   return s;
2815 }
2816 @
2817
2818 <<tension balance function tests>>=
2819 <<check relative error>>
2820
2821 double hooke(double x, void *pK)
2822 {
2823   assert(x >= 0);
2824   return *((double*)pK) * x;
2825 }
2826
2827 START_TEST(test_single_function)
2828 {
2829   double k=5, x=3, last_x=2.0, F;
2830   one_dim_fn_t *handlers[] = {&hooke};
2831   one_dim_fn_t *inverse_handlers[] = {NULL};
2832   void *data[] =             {&k};
2833   F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2834   fail_unless(F = k*x, NULL);
2835 }
2836 END_TEST
2837 @
2838
2839 We can also test balancing two springs with different spring constants.
2840 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2841 <<tension balance function tests>>=
2842 START_TEST(test_double_hooke)
2843 {
2844   double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2845   one_dim_fn_t *handlers[] = {&hooke, &hooke};
2846   one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2847   void *data[] =             {&k1,    &k2};
2848   F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2849   x1e = x*k2/(k1+k2);
2850   Fe = k1*x1e;
2851   x2e = x1e*k1/k2;
2852   //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2853   //CHECK_ERR(1e-6, x1e, xi[0]);
2854   //CHECK_ERR(1e-6, x2e, xi[1]);
2855   CHECK_ERR(1e-6, Fe, F);
2856 }
2857 END_TEST
2858 @
2859 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2860
2861 <<tension balance function test case defs>>=
2862 TCase *tc_tbfunc = tcase_create("tension balance function");
2863 @
2864
2865 <<tension balance function test case adds>>=
2866 tcase_add_test(tc_tbfunc, test_single_function);
2867 tcase_add_test(tc_tbfunc, test_double_hooke);
2868 suite_add_tcase(s, tc_tbfunc);
2869 @
2870
2871 \section{Lists}
2872 \label{sec.lists}
2873
2874 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2875 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2876 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2877
2878 <<list.h>>=
2879 #ifndef LIST_H
2880 #define LIST_H
2881 <<license comment>>
2882 <<list definitions>>
2883 <<list declarations>>
2884 #endif /* LIST_H */
2885 @
2886
2887 <<list declarations>>=
2888 <<head and tail declarations>>
2889 <<list length declaration>>
2890 <<push declaration>>
2891 <<pop declaration>>
2892 <<list destroy declaration>>
2893 <<list to array declaration>>
2894 <<move declaration>>
2895 <<sort declaration>>
2896 <<list index declaration>>
2897 <<list element declaration>>
2898 <<unique declaration>>
2899 <<list print declaration>>
2900 @
2901
2902 <<list functions>>=
2903 <<create and destroy node>>
2904 <<head and tail>>
2905 <<list length>>
2906 <<push>>
2907 <<pop>>
2908 <<list destroy>>
2909 <<list to array>>
2910 <<move>>
2911 <<sort>>
2912 <<list index>>
2913 <<list element>>
2914 <<unique>>
2915 <<list print>>
2916 @
2917
2918 <<list module makefile lines>>=
2919 NW_SAWSIM_MODS += list
2920 CHECK_BINS += check_list
2921 check_list_MODS =
2922 @
2923
2924 <<list definitions>>=
2925 typedef struct list_struct {
2926   struct list_struct *next;
2927   struct list_struct *prev;
2928   void *d; /* data */
2929 } list_t;
2930
2931 <<comparison function typedef>>
2932 @
2933
2934 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2935 <<head and tail declarations>>=
2936 list_t *head(list_t *list);
2937 list_t *tail(list_t *list);
2938 @
2939 <<head and tail>>=
2940 list_t *head(list_t *list)
2941 {
2942   if (list == NULL)
2943     return list;
2944   while (list->prev != NULL) {
2945     list = list->prev;
2946   }
2947   return list;
2948 }
2949
2950 list_t *tail(list_t *list)
2951 {
2952   if (list == NULL)
2953     return list;
2954   while (list->next != NULL) {
2955     list = list->next;
2956   }
2957   return list;
2958 }
2959 @
2960
2961 <<list length declaration>>=
2962 int list_length(list_t *list);
2963 @
2964 <<list length>>=
2965 int list_length(list_t *list)
2966 {
2967   int i;
2968   if (list == NULL)
2969     return 0;
2970   list = head(list);
2971   i = 1;
2972   while (list->next != NULL) {
2973     list = list->next;
2974     i += 1;
2975   }
2976   return i;
2977 }
2978 @
2979
2980 [[push]] inserts a node relative to the current position in the list
2981 without changing the current position.
2982 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2983 so we set it to that of the pushed domain.
2984 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2985 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2986 <<push declaration>>=
2987 void push(list_t **pList, void *data, int below);
2988 @
2989 <<push>>=
2990 void push(list_t **pList, void *data, int below)
2991 {
2992   list_t *list, *new_node;
2993   assert(pList != NULL);
2994   list = *pList;
2995   new_node = create_node(data);
2996   if (list == NULL)
2997     *pList = new_node;
2998   else if (below==0) { /* insert above */
2999     if (list->prev != NULL)
3000       list->prev->next = new_node;
3001     new_node->prev = list->prev;
3002     list->prev = new_node;
3003     new_node->next = list;
3004   } else {         /* insert below */
3005     if (list->next != NULL)
3006       list->next->prev = new_node;
3007     new_node->next = list->next;
3008     list->next = new_node;
3009     new_node->prev = list;
3010   }
3011 }
3012 @
3013
3014 [[pop]] removes the current domain node, moving the current position
3015 to the node after it, unless that node is [[NULL]], in which case move
3016 the current position to the node before it.
3017 <<pop declaration>>=
3018 void *pop(list_t **pList);
3019 @
3020 <<pop>>=
3021 void *pop(list_t **pList)
3022 {
3023   list_t *list, *popped;
3024   void *data;
3025   assert(pList != NULL);
3026   list = *pList;
3027   assert(list != NULL); /* not an empty list */
3028   popped = list;
3029   /* bypass the popped node */
3030   if (list->prev != NULL)
3031      list->prev->next = list->next;
3032   if (list->next != NULL) {
3033      list->next->prev = list->prev;
3034      *pList = list->next;
3035   } else
3036      *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
3037   data = popped->d;
3038   destroy_node(popped);
3039   return data;
3040 }
3041 @
3042
3043 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
3044 If the cleanup function is [[NULL]], no cleanup function is called.
3045 The nodes are not popped in any particular order.
3046 <<list destroy declaration>>=
3047 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
3048 @
3049 <<list destroy>>=
3050 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
3051 {
3052   list_t *list;
3053   void *data;
3054   assert(pList != NULL);
3055   list = *pList;
3056   *pList = NULL;
3057   assert(list != NULL); /* not an empty list */
3058   while (list != NULL) {
3059     data = pop(&list);
3060     if (destroy != NULL)
3061       destroy(data);
3062   }
3063 }
3064 @
3065
3066 [[list_to_array]] converts a list to an array of pointers, preserving order.
3067 <<list to array declaration>>=
3068 void list_to_array(list_t *list, int *length, void ***array);
3069 @
3070 <<list to array>>=
3071 void list_to_array(list_t *list, int *pLength, void ***pArray)
3072 {
3073   void **array;
3074   int i,length;
3075   assert(list != NULL);
3076   assert(pLength != NULL);
3077   assert(pArray != NULL);
3078
3079   length = list_length(list);
3080   array = (void **)malloc(sizeof(void *)*length);
3081   assert(array != NULL);
3082   list = head(list);
3083   i=0;
3084   while (list != NULL) {
3085     array[i++] = list->d;
3086     list = list->next;
3087   }
3088   *pLength = length;
3089   *pArray = array;
3090 }
3091 @
3092
3093 [[move]] moves the current element along the list in either direction.
3094 <<move declaration>>=
3095 void move(list_t **pList, int downstream);
3096 @
3097 <<move>>=
3098 void move(list_t **pList, int downstream)
3099 {
3100   list_t *list, *flipper;
3101   void *data;
3102   assert(pList != NULL);
3103   list = *pList;          /* a>B>c>d */
3104   assert(list != NULL); /* not an empty list */
3105   if (downstream == 0)
3106     flipper = list->prev; /* flipper is A */
3107   else
3108     flipper = list->next; /* flipper is C */
3109   /* ensure there is somewhere to go in the direction we're heading */
3110   assert(flipper != NULL);
3111   /* remove the current element */
3112   data = pop(&list);      /* data=B, a>C>d */
3113   /* and put it back in in the correct spot */
3114   list = flipper;
3115   if (downstream == 0) {
3116     push(&list, data, 0); /* b>A>c>d */
3117     list = list->prev;    /* B>a>c>d   */
3118   } else {
3119     push(&list, data, 1); /* a>C>b>d */
3120     list = list->next;    /* a>c>B>d */
3121   }
3122   *pList = list;
3123 }
3124 @
3125
3126 [[sort]] sorts a list from smallest to largest according to a user
3127 specified comparison.
3128 <<comparison function typedef>>=
3129 typedef int (comparison_fn_t)(void *A, void *B);
3130 @
3131 For example
3132 <<int comparison function>>=
3133 int int_comparison(void *A, void *B)
3134 {
3135   int a,b;
3136   a = *((int *)A);
3137   b = *((int *)B);
3138   if (a > b) return 1;
3139   else if (a == b) return 0;
3140   else return -1;
3141 }
3142
3143 @
3144 <<sort declaration>>=
3145 void sort(list_t **pList, comparison_fn_t *comp);
3146 @
3147 Here we use the bubble sort algorithm.
3148 <<sort>>=
3149 void sort(list_t **pList, comparison_fn_t *comp)
3150 {
3151   list_t *list;
3152   int stable=0;
3153   assert(pList != NULL);
3154   list = *pList;
3155   assert(list != NULL); /* not an empty list */
3156   while (stable == 0) {
3157     stable = 1;
3158     list = head(list);
3159     while (list->next != NULL) {
3160       if ((*comp)(list->d, list->next->d) > 0) {
3161         stable = 0;
3162         move(&list, 1 /* downstream */);
3163       } else
3164         list = list->next;
3165     }
3166   }
3167   *pList = list;
3168 }
3169
3170 @
3171
3172 [[list_index]] finds the location of [[data]] in [[list]].  Returns
3173 [[-1]] if [[data]] is not in [[list]].
3174 <<list index declaration>>=
3175 int list_index(list_t *list, void *data);
3176 @
3177
3178 <<list index>>=
3179 int list_index(list_t *list, void *data)
3180 {
3181   int i=0;
3182   list = head(list);
3183   while (list != NULL) {
3184     if (list->d == data) return i;
3185     list = list->next;
3186     i++;
3187   }
3188   return -1;
3189 }
3190
3191 @
3192
3193 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3194 <<list element declaration>>=
3195 void *list_element(list_t *list, int i);
3196 @
3197 <<list element>>=
3198 void *list_element(list_t *list, int i)
3199 {
3200   int j=0;
3201   list = head(list);
3202   while (list != NULL) {
3203     if (i == j) return list->d;
3204     list = list->next;
3205     j++;
3206   }
3207   assert(0==1);
3208 }
3209
3210 @
3211
3212
3213 [[unique]] removes duplicates from a sorted list according to a user
3214 specified comparison.  Don't do this unless you have other pointers
3215 any dynamically allocated data in your list, because [[unique]] just
3216 drops any non-unique data without freeing it.
3217 <<unique declaration>>=
3218 void unique(list_t **pList, comparison_fn_t *comp);
3219 @
3220 <<unique>>=
3221 void unique(list_t **pList, comparison_fn_t *comp)
3222 {
3223   list_t *list;
3224   assert(pList != NULL);
3225   list = *pList;
3226   assert(list != NULL); /* not an empty list */
3227   list = head(list);
3228   while (list->next != NULL) {
3229     if ((*comp)(list->d, list->next->d) == 0) {
3230       pop(&list);
3231     } else
3232       list = list->next;
3233   }
3234   *pList = list;
3235 }
3236 @
3237
3238 [[list_print]] prints a list to a [[FILE *]] stream.
3239 <<list print declaration>>=
3240 void list_print(FILE *file, list_t *list, char *list_name);
3241 @
3242 <<list print>>=
3243 void list_print(FILE *file, list_t *list, char *list_name)
3244 {
3245   fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3246   list = head(list);
3247   while (list != NULL) {
3248     fprintf(file, " %p", list->d);
3249     list = list->next;
3250   }
3251   fprintf(file, "\n");
3252 }
3253 @
3254
3255 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3256 <<create and destroy node>>=
3257 list_t *create_node(void *data)
3258 {
3259   list_t *ret = (list_t *)malloc(sizeof(list_t));
3260   assert(ret != NULL);
3261   ret->prev = NULL;
3262   ret->next = NULL;
3263   ret->d = data;
3264   return ret;
3265 }
3266
3267 void destroy_node(list_t *node)
3268 {
3269   if (node)
3270     free(node);
3271 }
3272 @
3273 The user must free the data pointed to by the node on their own.
3274
3275 \subsection{List implementation}
3276
3277 <<list.c>>=
3278 <<license comment>>
3279 <<list includes>>
3280 #include "list.h"
3281 <<list functions>>
3282 @
3283
3284 <<list includes>>=
3285 #include <assert.h> /* assert()                                */
3286 #include <stdlib.h> /* malloc(), free(), rand()                */
3287 #include <stdio.h>  /* fprintf(), stdout, FILE                 */
3288 #include "global.h" /* destroy_data_func_t */
3289 @
3290
3291 \subsection{List unit tests}
3292
3293 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3294 <<check-list.c>>=
3295 <<license comment>>
3296 <<list check includes>>
3297 <<list test suite>>
3298 <<main check program>>
3299 @
3300
3301 <<list check includes>>=
3302 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3303 #include <stdio.h>  /* FILE                          */
3304 <<check includes>>
3305 #include "global.h"
3306 #include "list.h"
3307 @
3308
3309 <<list test suite>>=
3310 <<push tests>>
3311 <<pop tests>>
3312 <<list suite definition>>
3313 @
3314
3315 <<list suite definition>>=
3316 Suite *test_suite (void)
3317 {
3318   Suite *s = suite_create ("list");
3319   <<push test case defs>>
3320   <<pop test case defs>>
3321
3322   <<push test case adds>>
3323   <<pop test case adds>>
3324   return s;
3325 }
3326 @
3327
3328 <<push tests>>=
3329 START_TEST(test_push)
3330 {
3331   list_t *list=NULL;
3332   int i, p, e, n=3;
3333   for (i=0; i<n; i++)
3334     push(&list, (void *)i, 1);
3335   fail_unless(list == head(list), NULL);
3336   fail_unless( (int)list->d == 0 );
3337   for (i=0; i<n; i++) {
3338     /* we expect 0, n-1, n-2, ..., 1 */
3339     if (i == 0) e = 0;
3340     else        e = n-i;
3341     fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3342   }
3343 }
3344 END_TEST
3345 @
3346
3347 <<push test case defs>>=
3348 TCase *tc_push = tcase_create("push");
3349 @
3350
3351 <<push test case adds>>=
3352 tcase_add_test(tc_push, test_push);
3353 suite_add_tcase(s, tc_push);
3354 @
3355
3356 <<pop tests>>=
3357 @
3358 <<pop test case defs>>=
3359 @
3360 <<pop test case adds>>=
3361 @
3362
3363 \section{Function string evaluation}
3364
3365 For the saddle-point approximated Kramers' model (Section \ref{sec.kramers}) we need the ability to evaluate user-supplied functions ($E(x)$, $x_{ts}(F)$, \ldots).
3366 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3367 The solution is to run a scripting language as a subprocess accessed via pipes.
3368 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3369
3370 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3371 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3372 Most of this is also POSIX-specific (as far as I know), so you'll have to do some legwork to port it to non-POSIX systems.
3373 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3374
3375 If you feel the benefits do \emph{not} outweigh the costs, or you are on a non-POSIX system (e.g.\ MS Windows without Cygwin), you should probably hardcode your functions in \lang.
3376 Then you can either statically or dynamically link to those hardcoded functions.
3377 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3378
3379 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3380 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3381
3382 <<string-eval.h>>=
3383 #ifndef STRING_EVAL_H
3384 #define STRING_EVAL_H
3385 <<license comment>>
3386 <<string eval setup declaration>>
3387 <<string eval function declaration>>
3388 <<string eval teardown declaration>>
3389 #endif /* STRING_EVAL_H */
3390 @
3391
3392 <<string eval module makefile lines>>=
3393 NW_SAWSIM_MODS += string_eval
3394 CHECK_BINS += check_string_eval
3395 check_string_eval_MODS =
3396 @
3397
3398 For an introduction to POSIX process control, see\\
3399  \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3400  \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3401  and of course, the relavent [[man]] pages.
3402
3403 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3404 [[execvp]] replaces the calling process' program with a new program.
3405 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3406 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3407  but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3408 The new program needs command line arguments, just like it would if you were running it from a shell.
3409 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3410 with the final array entry being a [[NULL]] pointer.
3411
3412 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3413 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3414 <<string eval subprocess definitions>>=
3415 #define SUBPROCESS "python"
3416 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3417 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3418 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3419 @
3420 The [[i]] option lets Python know that it should run in interactive mode.
3421 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3422 In interactive mode, python acts on every instruction as soon as it is recieved.
3423 The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply turns off the default prompting ([[ps1]] and [[ps2]]), since that is mostly for human users, and our program doesn't need it.
3424 %The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply [[pass]]es so that the silly Python header information is not printed.  We leave the prompts in, because we scan for them to determine when the output has completed.
3425
3426 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3427 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3428 [[fork]] returns two copies of the same program, executing the original code.
3429 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3430 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3431
3432 We communicate with the child (Python) process using \emph{pipes}, with one process writing data into one end of the pipe, and the other process reading the data out of the other end.
3433 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3434 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3435 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3436 <<string eval pipe definitions>>=
3437 #define PIPE_READ  0   /* the end you read from */
3438 #define PIPE_WRITE 1   /* the end you write to */
3439
3440 #define STDIN      0   /* initial index of stdin pair  */
3441 #define STDOUT     2   /* initial index of stdout pair */
3442
3443 @
3444 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3445
3446 As a finishing touch, we can promote the POSIX file descriptors ([[read]]/[[write]] interface) into the more familiar [[stdio.h]] \emph{streams} ([[fprintf]]/[[fgetc]] interface) using [[fdopen]], which creates a stream from an open file descriptor.
3447
3448 <<string eval setup declaration>>=
3449 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3450 @
3451 <<string eval setup definition>>=
3452 void string_eval_setup(FILE **pIn, FILE **pOut)
3453 {
3454   pid_t pid;
3455   int pfd[4]; /* file descriptors for stdin and stdout */
3456   int rc;
3457   assert(pipe(pfd+STDIN)  != -1); /* stdin pair  (in, out) */
3458   assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3459
3460   pid = fork(); /* split process into two copies */
3461
3462   if (pid == -1) { /* fork error */
3463     perror("fork error");
3464     exit(1);
3465   } else if (pid == 0) { /* child */
3466     close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input   */
3467     close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3468     dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin  (closes original stdin) */
3469     dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3470     execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3471     perror("exec error"); /* exec shouldn't return */
3472     _exit(1);
3473   } else { /* parent */
3474     close(pfd[STDIN+PIPE_READ]);   /* close stdin pipe output */
3475     close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3476     *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3477     if ( *pIn == NULL ) {
3478       perror("fdopen (in)");
3479       exit(1);
3480     }
3481     *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3482     if ( *pOut == NULL ) {
3483       perror("fdopen (out)");
3484       exit(1);
3485     }
3486   }
3487 }
3488 @
3489
3490 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3491 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3492 <<string eval function declaration>>=
3493 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3494 @
3495 <<string eval function definition>>=
3496 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3497 {
3498   int rc;
3499   rc = fprintf(in, "%s", input);
3500   assert(rc == strlen(input));
3501   fflush(in);
3502   fflush(out);
3503   alarm(1); /* set a one second timeout on the read */
3504   assert( fgets(output, buflen, out) != NULL );
3505   alarm(0); /* cancel the timeout */
3506   //fprintf(stderr, "eval: %s ----> %s", input, output);
3507 }
3508 @
3509 The [[alarm]] calls set and clear a timeout on the returned output.
3510 If the timeout expires, the process would get a [[SIGALRM]], but it doesn't have a [[SIGALRM]] handler, so it gets a [[SIGKILL]] and dies.
3511 This protects against invalid input for which a line of output is not printed to [[stdout]].
3512 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3513 If you are getting strange results, check your python code seperately. TODO, better error handling.
3514
3515 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3516 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3517 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3518 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3519 <<string eval teardown declaration>>=
3520 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3521 @
3522
3523 <<string eval teardown definition>>=
3524 void string_eval_teardown(FILE **pIn, FILE **pOut)
3525 {
3526   pid_t pid=0;
3527   int stat_loc;
3528
3529   /* redirect Python's stderr */
3530   fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3531   fflush(*pIn);
3532
3533   /* close pipes */
3534   assert( fclose(*pIn) == 0 );
3535   *pIn = NULL;
3536   assert( fclose(*pOut) == 0 );
3537   *pOut = NULL;
3538
3539   /* wait for python to exit */
3540   while (pid <= 0) {
3541     pid = wait(&stat_loc);
3542     if (pid < 0) {
3543       perror("pid");
3544     }
3545   }
3546
3547   /*
3548   if (WIFEXITED(stat_loc)) {
3549     printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3550   } else if (WIFSIGNALED(stat_loc)) {
3551     printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3552   }
3553   */
3554 }
3555 @
3556 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3557
3558 \subsection{String evaluation implementation}
3559
3560 <<string-eval.c>>=
3561 <<license comment>>
3562 <<string eval includes>>
3563 #include "string_eval.h"
3564 <<string eval internal definitions>>
3565 <<string eval functions>>
3566 @
3567
3568 <<string eval includes>>=
3569 #include <assert.h> /* assert()                    */
3570 #include <stdlib.h> /* NULL                        */
3571 #include <stdio.h>  /* fprintf(), stdout, fdopen() */
3572 #include <string.h> /* strlen()                    */
3573 #include <sys/types.h> /* pid_t                    */
3574 #include <unistd.h> /* pipe(), fork(), execvp(), alarm()    */
3575 #include <wait.h>   /* wait()                      */
3576 @
3577
3578 <<string eval internal definitions>>=
3579 <<string eval subprocess definitions>>
3580 <<string eval pipe definitions>>
3581 @
3582
3583 <<string eval functions>>=
3584 <<string eval setup definition>>
3585 <<string eval function definition>>
3586 <<string eval teardown definition>>
3587 @
3588
3589 \subsection{String evaluation unit tests}
3590
3591 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3592 <<check-string-eval.c>>=
3593 <<license comment>>
3594 <<string eval check includes>>
3595 <<string comparison definition and globals>>
3596 <<string eval test suite>>
3597 <<main check program>>
3598 @
3599
3600 <<string eval check includes>>=
3601 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3602 #include <stdio.h>  /* printf()                      */
3603 #include <assert.h> /* assert()                      */
3604 #include <string.h> /* strlen()                      */
3605 #include <signal.h> /* SIGKILL                       */
3606 <<check includes>>
3607 #include "string_eval.h"
3608 @
3609
3610 <<string eval test suite>>=
3611 <<string eval tests>>
3612 <<string eval suite definition>>
3613 @
3614
3615 <<string eval suite definition>>=
3616 Suite *test_suite (void)
3617 {
3618   Suite *s = suite_create ("string eval");
3619   <<string eval test case defs>>
3620
3621   <<string eval test case adds>>
3622   return s;
3623 }
3624 @
3625
3626 <<string eval tests>>=
3627 START_TEST(test_setup_teardown)
3628 {
3629   FILE *in, *out;
3630   string_eval_setup(&in, &out);
3631   string_eval_teardown(&in, &out);
3632 }
3633 END_TEST
3634 START_TEST(test_invalid_command)
3635 {
3636   FILE *in, *out;
3637   char input[80], output[80]={};
3638   string_eval_setup(&in, &out);
3639   sprintf(input, "print undefined_name__traceback_should_be_printed_to_stderr\n");
3640   string_eval(in, out, input, 80, output);
3641   string_eval_teardown(&in, &out);
3642 }
3643 END_TEST
3644 START_TEST(test_simple_eval)
3645 {
3646   FILE *in, *out;
3647   char input[80], output[80]={};
3648   string_eval_setup(&in, &out);
3649   sprintf(input, "print 3+4*5\n");
3650   string_eval(in, out, input, 80, output);
3651   fail_unless(STRMATCH(output,"23\n"), NULL);
3652   string_eval_teardown(&in, &out);
3653 }
3654 END_TEST
3655 START_TEST(test_multiple_evals)
3656 {
3657   FILE *in, *out;
3658   char input[80], output[80]={};
3659   string_eval_setup(&in, &out);
3660   sprintf(input, "print 3+4*5\n");
3661   string_eval(in, out, input, 80, output);
3662   fail_unless(STRMATCH(output,"23\n"), NULL);
3663   sprintf(input, "print (3**2 + 4**2)**0.5\n");
3664   string_eval(in, out, input, 80, output);
3665   fail_unless(STRMATCH(output,"5.0\n"), NULL);
3666   string_eval_teardown(&in, &out);
3667 }
3668 END_TEST
3669 START_TEST(test_eval_with_state)
3670 {
3671   FILE *in, *out;
3672   char input[80], output[80]={};
3673   string_eval_setup(&in, &out);
3674   sprintf(input, "print 3+4*5\n");
3675   fprintf(in, "A = 3\n");
3676   sprintf(input, "print A*3\n");
3677   string_eval(in, out, input, 80, output);
3678   fail_unless(STRMATCH(output,"9\n"), NULL);
3679   string_eval_teardown(&in, &out);
3680 }
3681 END_TEST
3682 @
3683 <<string eval test case defs>>=
3684 TCase *tc_string_eval = tcase_create("string_eval");
3685 @
3686 <<string eval test case adds>>=
3687 tcase_add_test(tc_string_eval, test_setup_teardown);
3688 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3689 tcase_add_test(tc_string_eval, test_simple_eval);
3690 tcase_add_test(tc_string_eval, test_multiple_evals);
3691 tcase_add_test(tc_string_eval, test_eval_with_state);
3692 suite_add_tcase(s, tc_string_eval);
3693 @
3694
3695
3696 \section{Accelerating function evaluation}
3697
3698 My first version-0.3 code was running very slowly.
3699 With the options suggested in the help
3700 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3701 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3702 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3703 That's only 15 calls per solution, so the search algorithm seems reasonable.
3704 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3705
3706 <<accel-k.h>>=
3707 #ifndef ACCEL_K_H
3708 #define ACCEL_K_H
3709 <<license comment>>
3710 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3711 void free_accels();
3712 #endif /* ACCEL_K_H */
3713 @
3714
3715 <<accel k module makefile lines>>=
3716 NW_SAWSIM_MODS += accel_k
3717 CHECK_BINS += check_accel_k
3718 check_accel_k_MODS = interp k_model parse string_eval tavl
3719 @
3720
3721 <<accel-k.c>>=
3722 <<license comment>>
3723 #include <assert.h>  /* assert()                */
3724 #include <stdlib.h>  /* realloc(), free(), NULL */
3725 #include "global.h"  /* environment_t           */
3726 #include "k_model.h" /* k_func_t                */
3727 #include "interp.h"  /* interp_*                */
3728 #include "accel_k.h"
3729
3730 typedef struct accel_k_struct {
3731   interp_table_t *itable;
3732   k_func_t *k;
3733   environment_t *env;
3734   void *params;
3735 } accel_k_t;
3736
3737 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3738 static int num_accels = 0;
3739 static accel_k_t *accels=NULL;
3740
3741 /* Wrap k in the standard f(x) acceleration form */
3742 static double k_wrap(double F, void *params)
3743 {
3744   accel_k_t *p = (accel_k_t *)params;
3745   return (*p->k)(F, p->env, p->params);
3746 }
3747
3748 static int k_tol(double FA, double kA, double FB, double kB)
3749 {
3750   assert(FB > FA);
3751   if (FB-FA > 1e-12) {
3752     //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3753     return 1; /* unnacceptable */
3754   } else {
3755     //printf("acceptable tol\n");
3756     return 0; /* acceptable */
3757   }
3758 }
3759
3760 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3761 {
3762   int i=num_accels;
3763   accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3764   assert(accels != NULL);
3765   accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3766   accels[i].k = k;
3767   accels[i].env = env;
3768   accels[i].params = params;
3769   return i;
3770 }
3771
3772 void free_accels()
3773 {
3774   int i;
3775   for (i=0; i<num_accels; i++)
3776     interp_table_free(accels[i].itable);
3777   num_accels=0;
3778   if (accels) free(accels);
3779   accels = NULL;
3780 }
3781
3782 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3783 {
3784   int i;
3785   for (i=0; i<num_accels; i++) {
3786     if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3787       /* We've already setup this function.
3788        * WARNING: we're assuming param and env are the same. */
3789       return interp_table_eval(accels[i].itable, accels+i, F);
3790     }
3791   }
3792   /* set up a new acceleration environment */
3793   i = add_accel_k(k, env, params);
3794   return interp_table_eval(accels[i].itable, accels+i, F);
3795 }
3796 @
3797
3798 \subsection{Unit tests}
3799
3800 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3801 <<check-accel-k.c>>=
3802 <<license comment>>
3803 <<accel k check includes>>
3804 <<check relative error>>
3805 <<model globals>>
3806 <<accel k test suite>>
3807 <<main check program>>
3808 @
3809
3810 <<accel k check includes>>=
3811 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
3812 #include <stdio.h>  /* sprintf()                  */
3813 #include <assert.h> /* assert()                   */
3814 #include <math.h>   /* exp()                      */
3815 #include <errno.h>  /* errno, ERANGE              */
3816 <<check includes>>
3817 #include "global.h"
3818 #include "k_model.h"
3819 #include "accel_k.h"
3820
3821 @
3822
3823 <<accel k test suite>>=
3824 <<accel k suite tests>>
3825 <<accel k suite definition>>
3826 @
3827
3828 <<accel k suite definition>>=
3829 Suite *test_suite (void)
3830 {
3831   Suite *s = suite_create ("accelerated k model");
3832   <<accel k test case defs>>
3833
3834   <<accel k test case adds>>
3835   return s;
3836 }
3837 @
3838
3839 <<accel k test case defs>>=
3840 TCase *tc_accel_k = tcase_create("accelerated k model");
3841 @
3842
3843 <<accel k test case adds>>=
3844 tcase_add_test(tc_accel_k, test_accel_k_accuracy);
3845 suite_add_tcase(s, tc_accel_k);
3846 @
3847
3848 <<accel k suite tests>>=
3849 START_TEST(test_accel_k_accuracy)
3850 {
3851   double k, k_accel, F, Fm=0.0, FM=300e-12, dF=0.5e-12;
3852   k_func_t *k_func = bell_k;
3853   char *param_strings[] = {"3.3e-4", "0.25e-9"};
3854   void *params = string_create_bell_param_t(param_strings);
3855   environment_t env = {};
3856   env.T = 300.0;
3857
3858   for(F=Fm; F <= FM; F+=dF) {
3859     k = k_func(F, &env, params);
3860     k_accel = accel_k(k_func, F, &env, params);
3861     CHECK_ERR(1e-8, k, k_accel);
3862   }
3863   destroy_bell_param_t(params);
3864 }
3865 END_TEST
3866
3867 @
3868
3869 \section{Tension models}
3870 \label{sec.tension_models}
3871
3872 TODO: tension model intro.
3873 See \url{http://web.mit.edu/course/3/3.11/www/pset03/Rec9.pdf} for a quick introduction to worm-like and freely-jointed chain models.
3874
3875 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3876 The interface is defined in [[tension_model.h]], the implementation in [[tension_model.c]], the unit testing in [[check_k_model.c]], and some other tidbits in [[tension_model_utils.c]].
3877
3878 <<tension-model.h>>=
3879 #ifndef TENSION_MODEL_H
3880 #define TENSION_MODEL_H
3881 <<license comment>>
3882 <<tension handler types>>
3883 <<constant tension model declarations>>
3884 <<hooke tension model declarations>>
3885 <<worm-like chain tension model declarations>>
3886 <<freely-jointed chain tension model declarations>>
3887 <<piston tension model declarations>>
3888 <<find tension definitions>>
3889 <<tension model global declarations>>
3890 #endif /* TENSION_MODEL_H */
3891 @
3892
3893 <<tension model module makefile lines>>=
3894 NW_SAWSIM_MODS += tension_model
3895 CHECK_BINS += check_tension_model
3896 check_tension_model_MODS = list tension_balance
3897 @
3898 <<tension model utils makefile lines>>=
3899 TENSION_MODEL_MODS = tension_model parse list tension_balance
3900 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3901         $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3902         $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3903 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3904         notangle -Rtension-model-utils.c $< > $@
3905 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3906         $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3907 clean_tension_model_utils :
3908         rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3909 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3910 case, the directories) as ‘order-only’ prerequisites.  The timestamp
3911 on these prerequisits does not effect whether the rules are executed.
3912 This is appropriate for directories, where we don't need to recompile
3913 after an unrelated has been added to the directory, but only when the
3914 source prerequisites change.  See the [[make]] documentation for more
3915 details
3916 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3917
3918 \subsection{Null}
3919 \label{sec.null_tension}
3920
3921 For unstretchable domains.
3922
3923 <<null tension model getopt>>=
3924 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3925 @
3926
3927 \subsection{Constant}
3928 \label{sec.const_tension}
3929
3930 An infinitely stretchable domain providing a constant tension.
3931 Obviously this cannot be inverted, so you can't put this domain in
3932 series with any other domains.  We include it mostly for testing
3933 purposes.
3934
3935 <<constant tension functions>>=
3936 <<constant tension function>>
3937 <<constant tension parameter create/destroy functions>>
3938 @
3939
3940 <<constant tension model declarations>>=
3941 <<constant tension function declaration>>
3942 <<constant tension parameter create/destroy function declarations>>
3943 <<constant tension model global declarations>>
3944