src/sawsim.nw

   1 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
   2 % we have our own preamble,
   3 % use `noweave -delay` to not print noweb's automatic one
   4 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
   5 \documentclass[letterpaper, 11pt]{article}
   6 \usepackage{noweb}
   7 \pagestyle{noweb}
   8 \noweboptions{smallcode}
   9
  10 \usepackage{url}    % needed for natbib for underscores in urls?
  11 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
  12 % breaklinks breaks long links
  13 % colorlinks colors link text instead of boxing it
  14 \usepackage{amsmath}  % for \text{}, \xrightarrow[sub]{super}
  15 \usepackage{amsthm}  % for theorem and proof environments
  16 \newtheorem{theorem}{Theorem}
  17
  18 \usepackage{doc}    % BibTeX
  19 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
  20 % super selects citations in superscript mode
  21 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
  22
  23 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
  24 \usepackage[utf8x]{inputenc}  % allow UTF-8 characters
  25
  26 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
  27 %\bibliographystyle{plain} % pick the bibliography style, includes dates
  28 \bibliographystyle{plainnat}
  29 \defcitealias{sw:noweb}{{\tt noweb}}
  30 \defcitealias{galassi05}{GSL}
  31 \defcitealias{sw:check}{{\tt check}}
  32 \defcitealias{sw:python}{Python}
  33
  34 \topmargin -1.0in
  35 \headheight 0.5in
  36 \headsep 0.5in
  37 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
  38 \oddsidemargin -0.5in
  39 \textwidth 7.5in
  40 \setlength{\parindent}{0pt}
  41 \setlength{\parskip}{5pt}
  42
  43 % For some reason, the \maketitle wants to go on it's own page
  44 % so we'll just hardcode the following in our first page.
  45 %\title{Sawsim: a sawtooth protein unfolding simulator}
  46 %\author{W.~Trevor King}
  47 %\date{\today}
  48
  49 \newcommand{\prog}{[[sawsim]]}
  50 \newcommand{\lang}{[[C]]}
  51 \newcommand{\p}[3]{\left#1 #2 \right#3}   % e.g. \p({complicated expr})
  52 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}}   % ordinal th
  53 \newcommand{\U}[1]{\ensuremath{\text{ #1}}}      % units: $1\U{m/s}$
  54
  55 % single spaced lists, from Alvin Alexander
  56 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
  57 \newenvironment{packed_item}{
  58 \begin{itemize}
  59   \setlength{\itemsep}{1pt}
  60   \setlength{\parskip}{0pt}
  61   \setlength{\parsep}{0pt}
  62 }{\end{itemize}}
  63
  64 \begin{document}
  65 \nwfilename{sawsim.nw}
  66 \nwbegindocs{0}
  67
  68 %\maketitle
  69 \begin{centering}
  70   Sawsim: a sawtooth protein unfolding simulator \\
  71   W.~Trevor King \\
  72   \today \\
  73 \end{centering}
  74 \bigskip
  75 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  76
  77 \section{Introduction}
  78
  79 The unfolding of multi-globular protein strings is a stochastic
  80 process.  In the \prog\ program, we use Monte Carlo methods to
  81 simulate this unfolding through an explicit time-stepping approach.
  82 We develop a program in \lang\ to simulate probability of unfolding
  83 under a constant extension velocity of a chain of globular domains.
  84
  85 In order to extract physical meaning from many mechanical unfolding
  86 simulations, it is necessary to compare the experimental results
  87 against the results of simulations using different models and
  88 parameters.  The model and parameters that best match the experimental
  89 data provide the ‘best’ model for that experiment.
  90
  91 Previous mechanical unfolding studies have rolled their own
  92 simulations for modelling their experiment, and there has been little
  93 exchange and standardization of these simulations between groups.
  94 The \prog\ program attempts to fill this gap, providing a flexible,
  95 extensible, \emph{community} program, to make it easier for groups
  96 to share the burden and increase the transparency of data analysis.
  97
  98 ‘Sawsim’ is blend of ‘sawtooth simulation’.
  99
 100 \subsection{Related work}
 101
 102 TODO References
 103
 104 \subsection{About this document}
 105
 106 This document was written in \citetalias{sw:noweb} with code blocks in
 107 \lang\ and comment blocks in \LaTeX.
 108
 109 \subsection{System Requirements and Dependencies}
 110
 111
 112 \subsection{Change Log}
 113
 114 Version 0.0 used only the unfolded domain WLC to determine the tension
 115 and had a fixed timestep $dt$.
 116
 117 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
 118 probability for a given domain was constant.
 119
 120 Version 0.2 added Kramers' model unfolding rate calculations, and a
 121 switch to select Bell's or Kramers' model.  This induced a major
 122 rewrite, introducing the tension group abstraction, and a split to
 123 multiple \lang\ source files.  Also added a unit-test suites for
 124 sanity checking, and switched to SI units throughout.
 125
 126 Version 0.3 added integrated Kramers' model unfolding rate
 127 calculations.  Also replaced some of my hand-coded routines with GNU
 128 Scientific Library \citepalias{galassi05} calls.
 129
 130 Version 0.4 added Python evaluation for the saddle-point Kramers'
 131 model unfolding rate calculations, so that the model functions could
 132 be controlled from the command line.  Also added interpolating lookup
 133 tables to accelerate slow unfolding rate model evaluation, and command
 134 line control of tension groups.
 135
 136 Version 0.5 added the stiffness environmental parameter and it's
 137 associated unfolding models.
 138
 139 <<version definition>>=
 140 #define VERSION "0.5"
 141 @
 142
 143 \subsection{License}
 144
 145 <<license>>=
 146  sawsim - program for simulating single molecule mechanical unfolding.
 147  Copyright (C) 2008, William Trevor King
 148
 149  This program is free software; you can redistribute it and/or
 150  modify it under the terms of the GNU General Public License as
 151  published by the Free Software Foundation; either version 3 of the
 152  License, or (at your option) any later version.
 153
 154  This program is distributed in the hope that it will be useful, but
 155  WITHOUT ANY WARRANTY; without even the implied warranty of
 156  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 157  See the GNU General Public License for more details.
 158
 159  You should have received a copy of the GNU General Public License
 160  along with this program; if not, write to the Free Software
 161  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 162  02111-1307, USA.
 163
 164  The author may be contacted at <wking@drexel.edu> on the Internet, or
 165  write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
 166  Philadelphia PA 19104, USA.
 167 @
 168
 169 <<license comment>>=
 170 /*
 171  <<license>>
 172  */
 173
 174 @
 175
 176 \section{Main}
 177
 178 The unfolding system is stored as a chain of \emph{domains}.  Domains
 179 can be folded, globular protein domains, unfolded protein linkers, AFM
 180 cantilevers, or other stretchable system components.
 181
 182 Each domain contributes to the total tension, which depends on the
 183 chain extension.  The particular model for the tension as a function
 184 of extension is handled generally with function pointers.  So far the
 185 following models have been implemented:
 186 \begin{packed_item}
 187  \item  Null (Appendix \ref{sec.null_tension}),
 188  \item  Constant (Appendix \ref{sec.const_tension}),
 189  \item  Hookean spring (Appendix \ref{sec.hooke}),
 190  \item  Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
 191  \item  Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
 192 \end{packed_item}
 193
 194 The tension model and parameters used for a particular domain depend
 195 on whether the domain is folded or unfolded.  The transition rate from
 196 the folded state to the unfolded state is also handled generally with
 197 function pointers, with implementations for
 198 \begin{packed_item}
 199  \item  Null (Appendix \ref{sec.null_k}),
 200  \item  Bell's model (Appendix \ref{sec.bell}),
 201  \item  Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
 202  \item  Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
 203  \item  Kramers' saddle point model (Appendix \ref{sec.kramers}).
 204 \end{packed_item}
 205
 206 The unfolding of the chain domains is modeled in two stages.  First
 207 the tension of the chain is calculated.  Then the tension (assumed to
 208 be constant over the length of the chain), is applied to each folded
 209 domain, and used to calculate the probability of that subdomain
 210 unfolding in the next timestep $dt$.  Then the time is stepped
 211 forward, and the process is repeated.
 212
 213 <<main program>>=
 214 int main(int argc, char **argv)
 215 {
 216   <<tension model getopt array definition>>
 217   <<k model getopt array definition>>
 218   list_t *domain_list=NULL;    /* variables initialized in get_options() */
 219   environment_t env={0};
 220   double P, dt_max, v, xstep;
 221   int num_folded=0, unfolding;
 222   double x=0, dt, F;           /* variables used in the simulation loop */
 223   get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS,
 224             tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
 225   setup();
 226   if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
 227   if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
 228   while (num_folded > 0) {
 229     F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
 230     if (xstep == 0)
 231       dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
 232     else
 233       dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, 0);
 234     unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
 235     if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
 236     if (xstep == 0) {
 237       x += v*dt;
 238     } else {
 239       if (dt == dt_max) { /* step completed */
 240         x += xstep;
 241         dt_max = xstep / v;
 242       } else { /* still working on this step */
 243         dt_max -= dt;
 244       }
 245     }
 246   }
 247   destroy_domain_list(domain_list);
 248   free_accels();
 249   return 0;
 250 }
 251 @ The meat of the simulation is bundled into the three functions
 252 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
 253 [[find_tension]] is discussed in Section \ref{sec.find_tension},
 254 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
 255 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
 256
 257 The stretched distance is extended in one of two modes depending on
 258 [[xstep]].  If [[xstep == 0]], then the pulling is continuous, at
 259 least within the limits of the inherent discretization of a time
 260 stepped approach.  At any rate, the timestep is limited by the maximum
 261 allowed timestep [[dt_max]] and the maximum allowed unfolding
 262 probability in a given step [[P]].  In the second mode [[xstep > 0]],
 263 and the end to end distance increases in discrete steps of that
 264 length.  The time between steps is chosen to maintain an average
 265 unfolding velocity [[v]].  This approach is less work to simulate than
 266 smooth pulling and also closer to the reality of many experiments, but
 267 it is practically impossible to treat analytically.  With both pulling
 268 styles implemented, the effects of the discretization can be easily
 269 calculated, bridging the gap between theory and experiment.
 270
 271 Environmental parameters important in determining reaction rates and
 272 tensions (e.g. temperature, pH) are stored in a single structure to
 273 facilitate extension to more complicated models in the future.
 274 <<environment definition>>=
 275 typedef struct environment_struct {
 276   double T;
 277   double stiffness;
 278 } environment_t;
 279 @
 280
 281 <<global.h>>=
 282 #ifndef GLOBAL_H
 283 #define GLOBAL_H
 284 <<environment definition>>
 285 <<one dimensional function definition>>
 286 <<create/destroy definitions>>
 287 <<model globals>>
 288 #endif /* GLOBAL_H */
 289 @
 290 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
 291 included multiple times.
 292
 293 \section{Simulation functions}
 294
 295 <<simulation functions>>=
 296 <<find tension>>
 297 <<determine dt>>
 298 <<happens>>
 299 <<does domain unfold>>
 300 <<randomly unfold domains>>
 301 @
 302
 303 \subsection{Tension}
 304 \label{sec.find_tension}
 305
 306 Because the stretched system may be made up of several parts (folded
 307 domains, unfolded domains, spring-like cantilever, \ldots), we will
 308 assign the domains to models and groups.  The models are used to
 309 determine a domain's tension (Hookean spring, WLC, \ldots).  Groups
 310 will represent collections of domains which the model should treat as
 311 a single entity.  A domain may belong to different groups or models
 312 depending on its state.  For example, a domain might be modeled as a
 313 freely-jointed chain when it is folded and as a worm-like chain when
 314 it is unfolded.  The domains are assumed to be commutative, so
 315 ordering is ignored.  The interactions between the groups are assumed
 316 to be linear, but the interactions between domains of the same group
 317 need not be.  This allows for non-linear group models such as th
 318 worm-like or freely-jointed chains.  Each model has a tension handler
 319 function, which gives the tension $F$ needed to stretch a given group
 320 of domains a distance $x$.
 321
 322 To understand the purpose of groups, consider a chain of proteins
 323 where two folded proteins $A$ and $B$ are modeled as WLCs with
 324 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
 325 modeled as WLCs with persistence length $p_u$.  The proteins are
 326 attached to a cantilever $E$ qof spring constant $κ$.  The string
 327 would get split into three lists:
 328 \begin{center}
 329   \begin{tabular}{llll}
 330     Model & Group &    List & Tension \\
 331     WLC   &     0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
 332     WLC   &     1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
 333     Hooke &     0 & $[E]$   & $F_\text{Hooke}(x, κ)$ \\
 334   \end{tabular}
 335 \end{center}
 336 Note that group numbers only matter \emph{within} model classes, since
 337 grouping domains with seperate models doesn't make sense.
 338
 339 <<find tension definitions>>=
 340 #define NUM_TENSION_MODELS 5
 341
 342 @
 343
 344 <<tension handler array global declaration>>=
 345 extern one_dim_fn_t *tension_handlers[];
 346 @
 347
 348 <<tension handler array global>>=
 349 one_dim_fn_t *tension_handlers[] = {
 350               NULL,
 351               &const_tension_handler,
 352               &hooke_handler,
 353               &wlc_handler,
 354               &fjc_handler,
 355               };
 356
 357 @
 358
 359 <<tension model global declarations>>=
 360 <<tension handler array global declaration>>
 361 @
 362
 363 <<tension handler types>>=
 364 typedef struct tension_handler_data_struct {
 365   /* int            num_domains; */
 366   /* domain_t      *domains;     */
 367   list_t        *group;
 368   environment_t *env;
 369   void          *persist;
 370 } tension_handler_data_t;
 371
 372 @
 373
 374 After sorting the chain into separate groups $G_i$, with tension
 375 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
 376 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
 377 \forall i,j$.  Note that there may be several groups within each model
 378 class.  [[find_tension]] is basically a wrapper to feed proper input
 379 into the [[tension_balance]] function.  [[unfolding]] should be set to
 380 0 if there were no unfoldings since the previous call to
 381 [[find_tension]].  While were messing with the tension handlers, we
 382 also grab a measure of the current linker stiffness for the
 383 environmental variable, which is needed by some of the unfolding rate
 384 models (e.g. the linker-stiffness corrected Bell model, Appendix
 385 \ref{sec.kbell}).
 386 <<find tension>>=
 387 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
 388 {
 389   static int num_active_groups=0;
 390   static one_dim_fn_t **per_group_handlers = NULL;
 391   static void **per_group_params = NULL;
 392   static double last_x;
 393   tension_handler_data_t *tdata;
 394   double F;
 395   int i;
 396
 397 #ifdef DEBUG
 398   fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
 399   list_print(stderr, domain_list, "domain list");
 400 #endif
 401
 402   if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
 403     /* free old statics */
 404     if (per_group_handlers != NULL)
 405       free(per_group_handlers);
 406     if (per_group_params != NULL) {
 407       for (i=0; i < num_active_groups; i++) {
 408         assert(per_group_params[i] != NULL);
 409         free(per_group_params[i]);
 410       }
 411       free(per_group_params);
 412     }
 413
 414     /* get new statics */
 415     get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
 416
 417     /* fill in the group handlers and params */
 418     for (i=0; i<num_active_groups; i++) {
 419       tdata = (tension_handler_data_t *) per_group_params[i];
 420 #ifdef DEBUG
 421       fprintf(stderr, "active group %d of %d", i, num_active_groups);
 422       list_print(stderr, tdata->group, " parameters");
 423 #endif
 424       tdata->env = env;
 425       /* tdata->persist continues unchanged */
 426     }
 427     last_x = -1.0;
 428   } /* else roll with the current statics */
 429
 430   F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
 431
 432   {
 433     double dx;
 434     double xlow;
 435     double Flow;
 436 #ifdef DEBUG
 437   fprintf(stderr, "finding linker stiffness at distance %g, tension %g\n", x, F);
 438 #endif
 439     if (last_x == -1.0)
 440       dx = x/200.0;
 441     else
 442       dx = (x-last_x);
 443     if (dx == 0) { /* e.g. if x == 0 */
 444       env->stiffness = 0;
 445     } else {
 446       xlow = x-dx;
 447       Flow = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, xlow);
 448       env->stiffness = (F-Flow)/dx;
 449     }
 450 #ifdef DEBUG
 451   fprintf(stderr, "  stiffness (%g-%g)/%g = %g\t(dx = %g = %g-%g)\n", F, Flow, dx, env->stiffness, dx, x, xlow);
 452 #endif
 453   }
 454
 455   last_x = x;
 456
 457   return F;
 458 }
 459 @ For the moment, we will restrict our group boundaries to a single
 460 dimension, so $\sum_i x_i = x$, our total extension (it is this
 461 restriction that causes the groups to interact linearly).  We'll also
 462 restrict our extensions to all be positive.  With these restrictions,
 463 the problem of balancing the tensions reduces to solving a system of
 464 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
 465 the number of active groups.  In general this can be fairly
 466 complicated, but without much loss of practicality we can restrict
 467 ourselves to strictly monotonically increasing, non-negative tension
 468 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
 469 simpler.  For example, it guarantees the existence of a unique, real
 470 solution for finite forces.  See Appendix \ref{app.tension_balance}
 471 for the tension balancing implementation.
 472
 473 Each group must define a parameter structure if the tension model is
 474 parametrized,
 475  a function to create the parameter structure,
 476  a function to destroy the parameter structure, and
 477  a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
 478 The tension-specific parameter structures are contained in the domain
 479 group field.  See the Section \ref{app.model_selection} for a
 480 discussion on the command line framework.  See the worm-like chain
 481 implementation for an example (Section \ref{sec.wlc}).
 482
 483 The major limitation of our approach is that all of our tension models
 484 are Markovian (memory-less), which is not really a problem for the
 485 chains (see \citet{evans99} for WLC thermalization rate calculations).
 486
 487 \subsection{Unfolding rate}
 488 \label{sec.unfolding_rate}
 489
 490 Each folded domain is modeled as unimolecular, first order reaction
 491 $$
 492   \text{Folded} \xrightarrow{k}  % in tex: X atop Y
 493   \text{Unfolded}
 494 $$
 495 With the rate constant $k$ defined by
 496 $$
 497   \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
 498 $$
 499 So the probability of a given protein unfolding in a short time $dt$
 500 is given by
 501 $$
 502   P = k \cdot dt.
 503 $$
 504
 505 <<does domain unfold>>=
 506 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
 507 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
 508   double k;
 509   k = accel_k(domain->k, F, env, domain->k_params);
 510   //(*domain->k)(F, env, domain->k_params);
 511   //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
 512   return happens(k*dt); /* dice roll for prob. k*dt event */
 513 }
 514 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
 515
 516 Only one domain can unfold in each timestep, because the timescale of
 517 a domain unfolding $dt_u$ is assumed to be less than the simulation
 518 timestep $dt$, so a domain will completely unfold in a single
 519 timestep.  We adapt our timesteps to keep the probability of a single
 520 domain unfolding low, so the probability of two domains unfolding in
 521 the same timestep is negligible.  This approach breaks down as the
 522 adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
 523 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
 524 shouldn't be a problem.  To reassure yourself, you can ask the
 525 simulator to print the smallest timestep that was used with TODO.
 526 <<randomly unfold domains>>=
 527 int random_unfoldings(list_t *domain_list, double tension,
 528                       double dt, environment_t *env,
 529                       int *num_folded_domains)
 530 { /* returns 1 if there was an unfolding and 0 otherwise */
 531   while (domain_list != NULL) {
 532     if (D(domain_list)->state == DS_FOLDED
 533         && domain_unfolds(tension, dt, env, D(domain_list))) {
 534       if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
 535         fprintf(stdout, "%g\n", tension);
 536       D(domain_list)->state = DS_UNFOLDED;
 537       (*num_folded_domains)--;
 538       return 1; /* our one domain has unfolded, stop looking for others */
 539     }
 540     domain_list = domain_list->next;
 541   }
 542   return 0;
 543 }
 544 @
 545
 546 \subsection{Adaptive timesteps}
 547 \label{sec.adaptive_dt}
 548
 549 We'd like to pick $dt$ so the probability of unfolding in the next
 550 timestep is small.  If we don't adapt our timestep, we risk breaking
 551 our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
 552 only one domain unfolds in a given timestep.  Because $F(x)$ is
 553 monotonically increasing, excessively large timesteps will lead to
 554 erroneously large unfolding forces.  The simulation would have been
 555 accurate for sufficiently small fixed timesteps, but adaptive
 556 timesteps will allow us to move through low tension regions in fewer
 557 steps, leading to a more efficient simulation.
 558
 559 The actual adaptive timestep implementation is not particularly
 560 interesting, since we are only required to reduce $dt$ to somewhere
 561 below a set threshold, so I've removed it to Appendix
 562 \ref{app.adaptive_dt}.
 563
 564 \section{Models}
 565
 566 TODO: model intro.
 567
 568 The models provide the physics of the simulation, but the simulation
 569 allows interchangeable models, and we are currently focusing on the
 570 simulation itself, so we remove the actual model implementation to the
 571 appendices.
 572
 573 The tension models are defined in Appendix \ref{sec.tension_models}
 574 and unfolding models are defined in Appendix \ref{sec.k_models}.
 575
 576 <<model globals>>=
 577 #define k_B   1.3806503e-23  /* J/K */
 578 @
 579
 580
 581 \section{Command line interface}
 582
 583 <<get option functions>>=
 584 <<help>>
 585 <<index tension model>>
 586 <<index k model>>
 587 <<generate domain>>
 588 <<get options>>
 589 @
 590
 591 \subsection{Model selection}
 592 \label{app.model_selection}
 593
 594 The main difficulty with the command line interface in \prog\ is
 595 developing an intuitive method for accessing the possibly large number
 596 of available models.  We'll treat this generally by defining an array
 597 of available models, containing their important parameters.
 598
 599 <<get options definitions>>=
 600 <<model getopt definitions>>
 601 @
 602
 603 <<create/destroy definitions>>=
 604 typedef void *create_data_func_t(char **param_strings);
 605 typedef void destroy_data_func_t(void *param_struct);
 606 @
 607
 608 <<model getopt definitions>>=
 609 <<tension model getopt definitions>>
 610 <<k model getopt definitions>>
 611 @
 612
 613
 614 \subsubsection{Tension}
 615
 616 To access the various tension models from the command line, we define
 617 a structure that contains all the neccessary information about the
 618 model.
 619 <<tension model getopt definitions>>=
 620 typedef struct tension_model_getopt_struct {
 621   one_dim_fn_t      *handler;     /* fn implementing the model on a group */
 622   char                *name;        /* name string identifying the model    */
 623   char                *description; /* a brief description of the model     */
 624   int                  num_params;  /* number of extra params for the model */
 625   char               **param_descriptions; /* descriptions of extra params  */
 626   char                *params;      /* default values for the extra params  */
 627   create_data_func_t  *creator;     /* param-string -> model-data-struct fn */
 628   destroy_data_func_t *destructor;  /* fn to free the model data structure  */
 629 } tension_model_getopt_t;
 630 @
 631
 632 <<tension model getopt array definition>>=
 633 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
 634   <<null tension model getopt>>,
 635   <<constant tension model getopt>>,
 636   <<hooke tension model getopt>>,
 637   <<worm-like chain tension model getopt>>,
 638   <<freely-jointed chain tension model getopt>>
 639 };
 640 @
 641
 642 \subsubsection{Unfolding rate}
 643
 644 <<k model getopt definitions>>=
 645 #define NUM_K_MODELS 6
 646
 647 typedef struct k_model_getopt_struct {
 648   char *name;
 649   char *description;
 650   k_func_t *k;
 651   int num_params;
 652   char **param_descriptions;
 653   char *params;
 654   create_data_func_t *creator;
 655   destroy_data_func_t *destructor;
 656 } k_model_getopt_t;
 657 @
 658
 659 <<k model getopt array definition>>=
 660 k_model_getopt_t k_models[NUM_K_MODELS] = {
 661   <<null k model getopt>>,
 662   <<const k model getopt>>,
 663   <<bell k model getopt>>,
 664   <<kbell k model getopt>>,
 665   <<kramers k model getopt>>,
 666   <<kramers integ k model getopt>>
 667 };
 668 @
 669
 670 \subsection{help}
 671
 672 <<help>>=
 673 void help(char *prog_name, double P, double dt_max, double v, double xstep,
 674           environment_t *env,
 675           int n_tension_models, tension_model_getopt_t *tension_models,
 676           int folded_tension_model, int unfolded_tension_model,
 677           int n_k_models, k_model_getopt_t *k_models,
 678           int k_model)
 679 {
 680   int i, j;
 681   printf("usage: %s [options]\n", prog_name);
 682   printf("Version %s\n\n", VERSION);
 683   printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
 684   printf("Simulation options:\n");
 685   printf("-P P\tTarget probability for dt (currently %g)\n", P);
 686   printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
 687   printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
 688   printf("-x dx\tPulling step size (currently %g m)\n", xstep);
 689   printf("\tWhen dx = 0, the pulling is continuous\n");
 690   printf("\tWhen dx > 0, the pulling occurs in discrete steps\n");
 691   printf("Environmental options:\n");
 692   printf("-T T\tTemperature T (currently %g K)\n", env->T);
 693   printf("-C T\tYou can also set the temperature T in Celsius\n");
 694   printf("Model options:\n");
 695   printf("The domains exist in either the folded or unfolded state\n");
 696   printf("The following options change the default behavior in each state.\n");
 697   printf("-m model[,group]\tFolded tension model (currently %s)\n",
 698          tension_models[folded_tension_model].name);
 699   printf("-a args\tFolded tension model argument string (currently %s)\n",
 700          tension_models[folded_tension_model].params);
 701   printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
 702          tension_models[unfolded_tension_model].name);
 703   printf("-A args\tUnfolded tension model argument string (currently %s)\n",
 704          tension_models[unfolded_tension_model].params);
 705   printf("The following options change the unfolding rate.\n");
 706   printf("-k model\tTransition rate model (currently %s)\n",
 707          k_models[k_model].name);
 708   printf("-K args\tTransition rate model argument string (currently %s)\n",
 709          k_models[k_model].params);
 710   printf("Domain creation options:\n");
 711   printf("Once you've set up the appropriate default models, you need to add the domains\n");
 712   printf("-F n\tAdd n folded domains with the current models\n");
 713   printf("-U n\tAdd n unfolded domains with the current models\n");
 714   printf("Output mode options:\n");
 715   printf("There are two output modes.  In standard mode, only the unfolding\n");
 716   printf("events are printed.  For example:\n");
 717   printf("  #Unfolding Force (N)\n");
 718   printf("  123.456e-12\n");
 719   printf("  ...\n");
 720   printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
 721   printf("  #Distance (m)\tForce (N)\tStiffness (N/m)\n");
 722   printf("  0.001\t0.0023\n");
 723   printf("  ...\n");
 724   printf("-V\tChange output to verbose mode\n");
 725   printf("-h\tPrint this help and exit\n");
 726   printf("\n");
 727   printf("Tension models:\n");
 728   for (i=0; i<n_tension_models; i++) {
 729     printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
 730     for (j=0; j < tension_models[i].num_params ; j++)
 731       printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
 732     printf("\t\tdefaults: %s\n", tension_models[i].params);
 733   }
 734   printf("Unfolding rate models:\n");
 735   for (i=0; i<n_k_models; i++) {
 736     printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
 737     for (j=0; j < k_models[i].num_params ; j++)
 738       printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
 739     printf("\t\tdefaults: %s\n", k_models[i].params);
 740   }
 741   printf("Examples:\n");
 742   printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
 743   printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
 744   printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
 745   printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
 746 }
 747 @
 748
 749 \subsection{Get options}
 750
 751 <<get options>>=
 752 void get_options(int argc, char **argv,
 753                  double *pP, double *pDt_max, double *pV, double *pXstep,
 754                  environment_t *env,
 755                  int n_tension_models, tension_model_getopt_t *tension_models,
 756                  int n_k_models, k_model_getopt_t *k_models,
 757                  list_t **pDomain_list, int *pNum_folded)
 758 {
 759   char *prog_name = NULL;
 760   char c, options[] = "P:t:v:x:T:C:m:a:M:A:k:K:F:U:Vh";
 761   int ftension_model=0, utension_model=0, k_model=0;
 762   char *ftension_params=NULL, *utension_params=NULL;
 763   int i, n, ftension_group=0, utension_group=0;
 764   int num_strings;
 765   char **strings;
 766   extern char *optarg;
 767   extern int optind, optopt, opterr;
 768
 769   assert (argc > 0);
 770
 771 #ifdef DEBUG
 772   for (i=0; i<n_tension_models; i++) {
 773     fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
 774             i, tension_models[i].name, tension_models[i].handler);
 775     assert(tension_models[i].handler == tension_handlers[i]);
 776   }
 777 #endif
 778
 779   /* setup defaults */
 780   flags = FLAG_OUTPUT_UNFOLDING_FORCES;
 781   prog_name = argv[0];
 782   *pP = 1e-3;       /* % pop per s */
 783   *pDt_max = 0.001; /* s           */
 784   *pV = 1e-6;       /* m/s         */
 785   *pXstep = 0.0;    /* m           */
 786   env->T = 300.0;   /* K           */
 787
 788   while ((c=getopt(argc, argv, options)) != -1) {
 789     switch(c) {
 790     case 'P':  *pP      = atof(optarg);           break;
 791     case 't':  *pDt_max = atof(optarg);           break;
 792     case 'v':  *pV      = atof(optarg);           break;
 793     case 'x':  *pXstep  = atof(optarg);           break;
 794     case 'T':  env->T   = atof(optarg);           break;
 795     case 'C':  env->T   = atof(optarg)+273.15;    break;
 796     case 'm':
 797       parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
 798       assert(num_strings == 1 || num_strings == 2);
 799       if (num_strings == 1)
 800         ftension_group = 0;
 801       else
 802         ftension_group = atoi(strings[1]);
 803       ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
 804       ftension_params = tension_models[ftension_model].params;
 805       free_parsed_list(num_strings, strings);
 806       break;
 807     case 'a':
 808       ftension_params = optarg;
 809       break;
 810     case 'M':
 811       parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
 812       assert(num_strings == 1 || num_strings == 2);
 813       if (num_strings == 1)
 814         utension_group = 0;
 815       else
 816         utension_group = atoi(strings[1]);
 817       utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
 818       utension_params = tension_models[utension_model].params;
 819       free_parsed_list(num_strings, strings);
 820       break;
 821     case 'A':
 822       utension_params = optarg;
 823       break;
 824     case 'k':
 825       k_model = index_k_model(n_k_models, k_models, optarg);
 826       break;
 827     case 'K':
 828       k_models[k_model].params = optarg;
 829       break;
 830     case 'F':
 831       n = atoi(optarg);
 832       assert(n > 0);
 833       for (i=0; i<n; i++) {
 834         push(pDomain_list, generate_domain(DS_FOLDED,
 835                                            tension_models+ftension_model,
 836                                            ftension_group,
 837                                            ftension_params,
 838                                            tension_models+utension_model,
 839                                            utension_group,
 840                                            utension_params,
 841                                            k_models+k_model), 1);
 842       }
 843       *pNum_folded += n;
 844       break;
 845     case 'U':
 846       n = atoi(optarg);
 847       assert(n > 0);
 848       for (i=0; i<n; i++) {
 849         push(pDomain_list, generate_domain(DS_UNFOLDED,
 850                                            tension_models+ftension_model,
 851                                            ftension_group,
 852                                            ftension_params,
 853                                            tension_models+utension_model,
 854                                            utension_group,
 855                                            utension_params,
 856                                            k_models+k_model), 1);
 857       }
 858       break;
 859     case 'V':
 860       flags = FLAG_OUTPUT_FULL_CURVE;
 861       break;
 862     case '?':
 863       fprintf(stderr, "unrecognized option '%c'\n", optopt);
 864       /* fall through to default case */
 865     default:
 866       help(prog_name, *pP, *pDt_max, *pV, *pXstep, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
 867       exit(1);
 868     }
 869   }
 870   /* check the arguments */
 871   assert(*pDomain_list != NULL); /* no domains! */
 872   assert(*pP > 0.0); assert(*pP < 1.0);
 873   assert(*pDt_max > 0.0);
 874   assert(*pV > 0.0);
 875   assert(env->T > 0.0);
 876   return;
 877 }
 878 @
 879
 880 <<index tension model>>=
 881 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
 882 {
 883   int i;
 884   for (i=0; i<n_models; i++)
 885     if (STRMATCH(models[i].name, name))
 886       break;
 887   if (i == n_models) {
 888     fprintf(stderr, "Unrecognized tension model '%s'\n", name);
 889     exit(1);
 890   }
 891   return i;
 892 }
 893 @
 894
 895 <<index k model>>=
 896 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
 897 {
 898   int i;
 899   for (i=0; i<n_models; i++)
 900     if (STRMATCH(models[i].name, name))
 901       break;
 902   if (i == n_models) {
 903     fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
 904     exit(1);
 905   }
 906   return i;
 907 }
 908 @
 909
 910 <<initialize model definition>>=
 911 /* requires int num_param_args and char **param_args in the current scope
 912  * usage:
 913  *  INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
 914  * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
 915  */
 916 #define INIT_MODEL(role, model, param_string, param_pointer) \
 917   do {                                                       \
 918     parse_list_string(param_string, SEP, '{', '}',           \
 919                       &num_param_args, &param_args);         \
 920     if (num_param_args != model->num_params) {               \
 921       fprintf(stderr,                                        \
 922               "%s model %s expected %d params,\n",           \
 923               role, model->name, model->num_params);         \
 924       fprintf(stderr,                                        \
 925               "not the %d params in '%s'\n",                 \
 926               num_param_args, param_string);                 \
 927       assert(num_param_args == model->num_params);           \
 928     }                                                        \
 929     if (model->creator)                                      \
 930       param_pointer = (*model->creator)(param_args);         \
 931     else                                                     \
 932       param_pointer = NULL;                                  \
 933     free_parsed_list(num_param_args, param_args);            \
 934   } while (0);
 935 @
 936
 937 <<generate domain>>=
 938 void *generate_domain(enum domain_state_t state,
 939                       tension_model_getopt_t *folded_model,
 940                       int folded_group,
 941                       char *folded_param_string,
 942                       tension_model_getopt_t *unfolded_model,
 943                       int unfolded_group,
 944                       char *unfolded_param_string,
 945                       k_model_getopt_t *k_model)
 946 {
 947   void *folded_params, *unfolded_params, *k_params;
 948   int num_param_args; /* for INIT_MODEL() */
 949   char **param_args;  /* for INIT_MODEL() */
 950
 951 #ifdef DEBUG
 952   fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
 953   fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
 954           k_model->name, k_model->params,
 955           folded_model->name, folded_model->handler, folded_group, folded_param_string,
 956           unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
 957 #endif
 958   INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
 959   INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
 960   INIT_MODEL("k", k_model, k_model->params, k_params);
 961   return create_domain(state,
 962                        k_model->k, k_params, k_model->destructor,
 963                        folded_model->handler, folded_group, folded_params, folded_model->destructor,
 964                        unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
 965                        );
 966 }
 967 @
 968
 969 \phantomsection
 970 \appendix
 971 \addcontentsline{toc}{section}{Appendicies}
 972
 973 \section{sawsim.c details}
 974
 975 \subsection{Layout}
 976
 977 The general layout of our simulation code is:
 978 <<*>>=
 979 <<license comment>>
 980 <<includes>>
 981 <<definitions>>
 982 <<globals>>
 983 <<functions>>
 984 <<main program>>
 985 @
 986
 987 We include [[math.h]], so don't forget to link to the libm with `-lm'.
 988 <<includes>>=
 989 #include <assert.h> /* assert()                                */
 990 #include <stdlib.h> /* malloc(), free(), rand()                */
 991 #include <stdio.h>  /* fprintf(), stdout                       */
 992 #include <string.h> /* strlen, strtok()                        */
 993 #include <math.h>   /* exp(), M_PI, sqrt()                     */
 994 #include <sys/types.h> /* pid_t (returned by getpid())         */
 995 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
 996 #include "global.h"
 997 #include "list.h"
 998 #include "tension_balance.h"
 999 #include "k_model.h"
1000 #include "tension_model.h"
1001 #include "parse.h"
1002 #include "accel_k.h"
1003 @
1004
1005 <<definitions>>=
1006 <<version definition>>
1007 <<flag definitions>>
1008 <<max/min definitions>>
1009 <<string comparison definition and globals>>
1010 <<initialize model definition>>
1011 <<get options definitions>>
1012 <<domain definitions>>
1013 @
1014
1015 <<globals>>=
1016 <<flag globals>>
1017 <<model globals>>
1018 @
1019
1020 <<functions>>=
1021 <<create/destroy domain>>
1022 <<destroy domain list>>
1023 <<group functions>>
1024 <<simulation functions>>
1025 <<boilerplate functions>>
1026 @
1027
1028 <<boilerplate functions>>=
1029 <<setup>>
1030 <<get option functions>>
1031 @
1032
1033 <<setup>>=
1034 void setup(void)
1035 {
1036   srand(getpid()*time(NULL)); /* seed rand() */
1037 }
1038 @
1039
1040 <<flag definitions>>=
1041 /* in octal b/c of prefixed '0' */
1042 #define FLAG_OUTPUT_FULL_CURVE       01
1043 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
1044 @
1045
1046 <<flag globals>>=
1047 static unsigned long int flags = 0;
1048 @
1049
1050 \subsection{Utilities}
1051 \label{app.utils}
1052
1053 <<max/min definitions>>=
1054 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1055 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1056 @
1057
1058 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1059 <<string comparison definition and globals>>=
1060 // Check if two strings match, return 1 if they do
1061 static char *temp_string_A;
1062 static char *temp_string_B;
1063 #define STRMATCH(a,b)   (temp_string_A=a, temp_string_B=b, \
1064         strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1065         !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1066         /* +1 to also compare the '\0' */
1067 @
1068
1069 We also define a macro for our [[check]] unit testing
1070 <<check relative error>>=
1071 #define CHECK_ERR(max_err, expected, received)                               \
1072   do {                                                                       \
1073     fail_unless( (received-expected)/expected < max_err,                     \
1074                  "relative error %g >= %g in %s (Expected %g, received %g)", \
1075                  (received-expected)/expected, max_err, #received,           \
1076                  expected, received);                                        \
1077     fail_unless(-(received-expected)/expected < max_err,                     \
1078                  "relative error %g >= %g in %s (Expected %g, received %g)", \
1079                 -(received-expected)/expected, max_err, #received,           \
1080                  expected, received);                                        \
1081   } while(0)
1082 @
1083
1084 <<happens>>=
1085 int happens(double probability)
1086 {
1087   assert(probability >= 0.0); assert(probability <= 1.0);
1088   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*/
1089 }
1090 @
1091
1092
1093 \subsection{Adaptive timesteps}
1094 \label{app.adaptive_dt}
1095
1096 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1097 so basing the timestep on the the unfolding probability at the current tension
1098 is dangerous, and we need to search for a $dt$ for which
1099 $P(F(x+v*dt)) < P_\text{target}$.
1100 There are two cases to consider.
1101 In the most common, no domains have unfolded since the last step,
1102 and we expect the next step to be slightly shorter than the previous one.
1103 In the less common, domains did unfold in the last step,
1104 and we expect the next step to be considerably longer than the previous one.
1105 <<search dt>>=
1106 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1107                  list_t *domain_list,
1108                  environment_t *env, double x,
1109                  double target_prob, double max_dt, double v)
1110 { /* Returns the timestep dt in seconds for the current folded domain.
1111    * Takes a list of tension handlers, the list of domains,
1112    * a pointer env to the environmental data, a starting separation x in m,
1113    * a target_prob between 0 and 1,
1114    * max_dt in s, stretching velocity v in m/s.
1115    */
1116   double F, k, dtCur, dtU, dtUCur, dtL, dt;
1117
1118   /* get upper bound using the current position */
1119   F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1120   //printf("Start with x = %g (F = %g)\n", x, F);
1121   k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1122   //printf("x %g\tF %g\tk %g\n", x, F, k);
1123   dtU = target_prob / k;    /* P = k dt, dtU is an upper bound on dt */
1124   if (dtU > max_dt) {
1125     //printf("overshot max_dt\n");
1126     dtU = max_dt;
1127   }
1128   if (v == 0) /* discrete stepping, so force is temporatily constant */
1129     return dtU;
1130
1131   /* set a lower bound on dt too */
1132   dtL = 0.0;
1133
1134   /* The dt determined above may produce illegitimate forces or ks.
1135    * Reduce the upper bound until we have valid ks. */
1136   dt = dtU;
1137   F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1138   while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1139     dtU /= 2.0;
1140     dt = dtU;
1141     F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1142   }
1143   //printf("Try for dt = %g (F = %g)\n", dt, F);
1144   k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1145   /* returned k may be -1.0 */
1146   //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1147   while (k == -1.0) { /* reduce step until we hit a valid k */
1148     dtU /= 2.0;
1149     dt = dtU; /* hopefully, we can use the max dt, see if it works */
1150     F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1151     //printf("Try for dt = %g (F = %g)\n", dt, F);
1152     k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1153     //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1154   }
1155   assert(dtU > 1e-14);      /* timestep to valid k too small */
1156   dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1157   if (dtUCur >= dt)
1158     return dt; /* dtU is safe. */
1159
1160   /* dtU wasn't safe, lets see what would be. */
1161   while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1162     dt = (dtU + dtL) / 2.0;
1163     F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1164     //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1165     k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1166     dtCur = target_prob / k;
1167     //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1168     if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1169       dtL = dt;
1170     else if (dtCur < dt) {  /* unsafe timestep back from x+dt, but...    */
1171       dtU = dt;             /* ... stepping out only dtCur would be safe */
1172       dtUCur = dtCur;
1173     } else
1174       break; /* dtCur = dt */
1175   }
1176   return MAX(dtUCur, dtL);
1177 }
1178 @
1179
1180 To determine $dt$ for an array of potentially different folded domains,
1181 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1182 <<determine dt>>=
1183 <<search dt>>
1184 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1185                     environment_t *env, double x,
1186                     double target_prob, double dt_max, double v)
1187 { /* Returns the timestep dt in seconds.
1188    * Takes the list of folded domains, target_prob between 0 and 1,
1189    * F in N, and T in K. */
1190   double dt=dt_max, new_dt;
1191   assert(target_prob > 0.0); assert(target_prob < 1.0);
1192   assert(dt_max > 0.0);
1193
1194   while (domain_list != NULL) {
1195     if (D(domain_list)->state == DS_FOLDED) {
1196       new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1197       dt = MIN(dt, new_dt);
1198     }
1199     domain_list = domain_list->next;
1200   }
1201   return dt;
1202 }
1203 @
1204
1205 \subsection{Domain data}
1206
1207 Currently domains exist in two states, folded and unfolded, and the
1208 only allowed transitions are folded $\rightarrow$ unfolded.  Of
1209 course, it wouldn't be too complicated to extent this to a multi-state
1210 system, with an array containing the domains group for each possible
1211 state, and a matrix of transition-rate-calculating functions.
1212 However, at this point such generality seems unnecessary at this
1213 point.
1214 <<domain definitions>>=
1215 enum domain_state_t {DS_FOLDED,
1216                      DS_UNFOLDED
1217 };
1218
1219 typedef struct domain_struct {
1220   enum domain_state_t state;
1221   one_dim_fn_t       *folded_handler;
1222   int                   folded_group;
1223   one_dim_fn_t       *unfolded_handler;
1224   int                   unfolded_group;
1225   k_func_t *k;           /* function returning unfolding rate */
1226   void *folded_params;   /* pointer to folded parameters   */
1227   void *unfolded_params; /* pointer to unfolded parameters */
1228   void *k_params;        /* pointer to k parameters        */
1229   destroy_data_func_t *destroy_folded;
1230   destroy_data_func_t *destroy_unfolded;
1231   destroy_data_func_t *destroy_k;
1232 } domain_t;
1233
1234 /* get the domain data for the current list node */
1235 #define D(list) ((domain_t *)(list)->d)
1236 /* get the tension params for the current list node */
1237 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1238                      ? ((domain_t *)(list)->d)->folded_params   \
1239                      : ((domain_t *)(list)->d)->unfolded_params)
1240 @
1241 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1242 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1243 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1244 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1245 We store them with the domain data so that [[destroy_domain]] doesn't have to know which type of domain it's cleaning up after.
1246
1247 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1248 <<create/destroy domain>>=
1249 domain_t *create_domain(enum domain_state_t state,
1250                         k_func_t *k,
1251                         void *k_params,
1252                         destroy_data_func_t *destroy_k,
1253                         one_dim_fn_t *folded_handler,
1254                         int folded_group,
1255                         void *folded_params,
1256                         destroy_data_func_t *destroy_folded,
1257                         one_dim_fn_t *unfolded_handler,
1258                         int unfolded_group,
1259                         void *unfolded_params,
1260                         destroy_data_func_t *destroy_unfolded)
1261 {
1262   domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1263   assert(ret != NULL);
1264   if (state == DS_FOLDED) {
1265     assert(k != NULL);  /* the pointer points somewhere valid  */
1266     assert(*k != NULL); /* and there is something useful there */
1267   } else
1268     assert(state == DS_UNFOLDED);
1269   ret->state = state;
1270   ret->folded_handler = folded_handler;
1271   ret->folded_group = folded_group;
1272   ret->unfolded_handler = unfolded_handler;
1273   ret->unfolded_group = unfolded_group;
1274   ret->k = k;
1275   ret->folded_params = folded_params;
1276   ret->unfolded_params = unfolded_params;
1277   ret->k_params = k_params;
1278   ret->destroy_folded = destroy_folded;
1279   ret->destroy_unfolded = destroy_unfolded;
1280   ret->destroy_k = destroy_k;
1281 #ifdef DEBUG
1282   fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1283 #endif
1284   return ret;
1285 }
1286
1287 void destroy_domain(domain_t *domain)
1288 {
1289   if (domain) {
1290     //printf("domain %p & %p\n", *domain, domain);
1291     if (domain->destroy_folded)
1292       (*domain->destroy_folded)(domain->folded_params);
1293     if (domain->destroy_unfolded)
1294       (*domain->destroy_unfolded)(domain->unfolded_params);
1295     if (domain->destroy_k)
1296       (*domain->destroy_k)(domain->k_params);
1297     free(domain);
1298   }
1299 }
1300 @
1301
1302 <<destroy domain list>>=
1303 void destroy_domain_list(list_t *domain_list)
1304 {
1305   domain_list = head(domain_list);
1306   while (domain_list != NULL)
1307     destroy_domain((domain_t *) pop(&domain_list));
1308 }
1309 @
1310
1311 \subsection{Domain model and group handling}
1312
1313 <<group functions>>=
1314 <<get tension handler>>
1315 <<get group>>
1316 <<int comparison function>>
1317 <<find possible groups>>
1318 <<is group active>>
1319 <<get group list>>
1320 <<get active groups>>
1321 @
1322
1323 <<get tension handler>>=
1324 one_dim_fn_t *get_tension_handler(domain_t *domain)
1325 {
1326   if (domain->state == DS_FOLDED)
1327     return domain->folded_handler;
1328   else {
1329     if (domain->state != DS_UNFOLDED)
1330       fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1331     assert(domain->state == DS_UNFOLDED);
1332     return domain->unfolded_handler;
1333   }
1334 }
1335 @
1336
1337 <<get group>>=
1338 int get_group(domain_t *domain)
1339 {
1340   if (domain->state == DS_FOLDED)
1341     return domain->folded_group;
1342   else {
1343     if (domain->state != DS_UNFOLDED)
1344       fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1345     assert(domain->state == DS_UNFOLDED);
1346     return domain->unfolded_group;
1347   }
1348 }
1349 @
1350
1351 We already know all possible tension classes, since they are all
1352 hardcoded into \prog.  However, there may be any number of possible
1353 groups.  We define a function [[find_possible_groups]] to search for
1354 possible (not neccessarily active) groups.  It can search for
1355 subgroups of a single [[handler]], or by setting [[handler]] to
1356 [[NULL]], subgroups of \emph{any} handler.  It returns a list of group
1357 integers.
1358 <<find possible groups>>=
1359 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1360   list_t *ret = NULL;
1361   list = head(list);
1362   while (list != NULL) {
1363     if (handler == NULL || get_tension_handler(D(list)) == handler) {
1364       push(&ret, &D(list)->folded_group, 1);
1365       push(&ret, &D(list)->unfolded_group, 1);
1366     }
1367     list = list->next;
1368   }
1369
1370   if (ret == NULL)
1371     return ret; /* no groups with this handler, no processing needed */
1372
1373   /* sort the ret list, and remove duplicates */
1374   sort(&ret, &int_comparison);
1375   unique(&ret, &int_comparison);
1376 #ifdef DEBUG
1377   fprintf(stderr, "handler %p possible groups:", handler);
1378   list = head(ret);
1379   while (list != NULL) {
1380     fprintf(stderr, " %d", *((int *)list->d));
1381     list = list->next;
1382   }
1383   fprintf(stderr, "\n");
1384 #endif
1385   return ret;
1386 }
1387 @
1388
1389 Search a [[list]] of domains, and determine whether a particular model
1390 class and group number combination is active.
1391 <<is group active>>=
1392 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1393 {
1394   list = head(list);
1395   while (list != NULL) {
1396     if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1397       return 1;
1398     list = list->next;
1399   }
1400   return 0;
1401 }
1402 @
1403
1404 Search a [[list]] of domains, and return all domains matching a
1405 particular model class and group number.
1406 <<get group list>>=
1407 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1408 {
1409   list_t *ret = NULL;
1410   list = head(list);
1411   while (list != NULL) {
1412     if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1413       push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1414 #ifdef DEBUG
1415       fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
1416 #endif
1417     }
1418     list = list->next;
1419   }
1420   return ret;
1421 }
1422 @
1423 Because all the node data in lists returned by [[get_group_list]] is
1424 also in the main domain list, you shouldn't destroy the node data
1425 popped off when destroying the group lists.  It will all get cleaned
1426 up when the main domain list is destroyed.
1427
1428 Put all this together to scan through a list of domains and construct
1429 an array of all the active groups.
1430 <<get active groups>>=
1431 void get_active_groups(list_t *list,
1432                      int num_tension_handlers, one_dim_fn_t **tension_handlers,
1433                      int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1434 {
1435   void **active_groups=NULL;
1436   one_dim_fn_t *handler, **per_group_handlers=NULL;
1437   int i, num_active_groups, *group;
1438   list_t *possible_group_numbers,  *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1439
1440   /* get all the active groups in a list */
1441   for (i=0; i<num_tension_handlers; i++) {
1442     handler = tension_handlers[i];
1443     if (handler == NULL) continue; /* tensionless handler */
1444     possible_group_numbers = head(find_possible_groups(list, handler));
1445     while (possible_group_numbers != NULL) {
1446       group = (int *)pop(&possible_group_numbers);
1447       if (is_group_active(list, handler, *group) == 1) {
1448         active_group_list = get_group_list(list, handler, *group);
1449         push(&active_groups_list, active_group_list, 1);
1450         push(&per_group_handlers_list, handler, 1);
1451 #ifdef DEBUG
1452         fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1453         list_print(stderr, active_group_list, "active group");
1454 #endif
1455       }
1456     }
1457   }
1458
1459   /* allocate the array we'll be returning */
1460   num_active_groups = list_length(active_groups_list);
1461   active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1462   assert(active_groups != NULL);
1463   per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1464   assert(per_group_handlers != NULL);
1465
1466   /* populate the array we'll be returning */
1467   active_groups_list = head(active_groups_list);
1468   for (i=0; i<num_active_groups; i++) {
1469     handler = pop(&per_group_handlers_list);
1470     per_group_handlers[i] = handler;
1471
1472     active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1473     assert(active_groups[i] != NULL);
1474     active_group_list = pop(&active_groups_list);
1475     ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1476     ((tension_handler_data_t *)active_groups[i])->env = NULL;
1477     ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1478   }
1479   assert(active_groups_list == NULL);
1480   assert(per_group_handlers_list == NULL);
1481
1482   *pNum_active_groups = num_active_groups;
1483   *pPer_group_handlers = per_group_handlers;
1484   *pActive_groups = active_groups;
1485 }
1486 @
1487
1488
1489 \section{String parsing}
1490
1491 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1492 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1493 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1494 need to take care of parsing those parameters themselves.
1495 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1496
1497 <<parse.h>>=
1498 #ifndef PARSE
1499 #define PARSE
1500 <<license comment>>
1501 <<parse definitions>>
1502 <<parse declarations>>
1503 #endif /* PARSE */
1504 @
1505
1506 <<parse module makefile lines>>=
1507 NW_SAWSIM_MODS += parse
1508 CHECK_BINS += check_parse
1509 check_parse_MODS =
1510 @
1511
1512 <<parse definitions>>=
1513 #define SEP ',' /* argument separator character */
1514 @
1515
1516 <<parse declarations>>=
1517 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1518                        int *num, char ***string_array);
1519 extern void free_parsed_list(int num, char **string_array);
1520 @
1521
1522 [[parse_list_string]] allocates memory, don't forget to free it
1523 afterward with [[free_parsed_list]].  It does not alter the original.
1524
1525 The string may start off with a [[deeper]] character
1526 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1527 [[x,y]], where the pointer is one character in on the copied string.
1528 However, when we go to free the memory, we need a pointer to the
1529 beginning of the string.  In order to accommodate this for a string
1530 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1531 the first $N$ elements point to the separated arguments, and let the
1532 last element point to the start of the copied string regardless of
1533 braces.
1534 <<parse delimited list string functions>>=
1535 /* TODO, split out into parse.hc */
1536 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1537 {
1538   int i=0, depth = 0;
1539   while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1540     if (string[i] == deeper) {depth++;}
1541     else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1542     i++;
1543   }
1544   return i;
1545 }
1546
1547 void parse_list_string(char *string, char sep, char deeper, char shallower,
1548                        int *num, char ***string_array)
1549 {
1550   char *str=NULL, **ret=NULL;
1551   int i, j, n;
1552   if (string==NULL || strlen(string) == 0) {    /* handle the trivial cases */
1553     *num = 0;
1554     *string_array = NULL;
1555     return;
1556   }
1557   /* make a copy of the string, so we don't change the original */
1558   str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1559   assert(str != NULL);
1560   strcpy(str, string); /* we know str is long enough */
1561   /* count the number of regions, so we can allocate pointers to them */
1562   i=-1; n=0;
1563   do {
1564     n++; i++; /* move on to next argument */
1565     i += next_delim_index(str+i, sep, deeper, shallower);
1566     //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1567   } while (str[i] != '\0');
1568   ret = (char **)malloc(sizeof(char *)*(n+1));
1569   assert(ret != NULL);
1570   /* replace the separators with '\0' & assign pointers */
1571   ret[n] = str; /* point to the front of the copied string */
1572   j=0;
1573   ret[0] = str;
1574   for(i=1; i<n; i++) {
1575     j += next_delim_index(str+j, sep, deeper, shallower);
1576     str[j++] = '\0';
1577     ret[i] = str+j; /* point to the separated arguments */
1578   }
1579   /* strip off bounding braces, if any */
1580   for(i=0; i<n; i++) {
1581     if (ret[i][0] == deeper) {
1582       ret[i][0] = '\0';
1583       ret[i]++;
1584     }
1585     if (ret[i][strlen(ret[i])-1] == shallower)
1586       ret[i][strlen(ret[i])-1] = '\0';
1587   }
1588   *num = n;
1589   *string_array = ret;
1590 }
1591
1592 void free_parsed_list(int num, char **string_array)
1593 {
1594   if (string_array) {
1595     /* frees all items, since they were allocated as a single string*/
1596     free(string_array[num]);
1597     /* and free the array of pointers */
1598     free(string_array);
1599   }
1600 }
1601 @
1602
1603 \subsection{Parsing implementation}
1604
1605 <<parse.c>>=
1606 <<license comment>>
1607 <<parse includes>>
1608 #include "parse.h"
1609 <<parse delimited list string functions>>
1610 @
1611
1612 <<parse includes>>=
1613 #include <assert.h> /* assert()                */
1614 #include <stdlib.h> /* NULL                    */
1615 #include <stdio.h>  /* fprintf(), stdout       *//*!!*/
1616 #include <string.h> /* strlen()                */
1617 #include "parse.h"
1618 @
1619
1620 \subsection{Parsing unit tests}
1621
1622 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1623 <<check-parse.c>>=
1624 <<parse check includes>>
1625 <<string comparison definition and globals>>
1626 <<check relative error>>
1627 <<parse test suite>>
1628 <<main check program>>
1629 @
1630
1631 <<parse check includes>>=
1632 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1633 #include <stdio.h>  /* printf()                              */
1634 #include <assert.h> /* assert()                              */
1635 #include <string.h> /* strlen()                              */
1636 <<check includes>>
1637 #include "parse.h"
1638 @
1639
1640 <<parse test suite>>=
1641 <<parse list string tests>>
1642 <<parse suite definition>>
1643 @
1644
1645 <<parse suite definition>>=
1646 Suite *test_suite (void)
1647 {
1648   Suite *s = suite_create ("k model");
1649   <<parse list string test case defs>>
1650
1651   <<parse list string test case add>>
1652   return s;
1653 }
1654 @
1655
1656 <<parse list string tests>>=
1657
1658 /*
1659 START_TEST(test_next_delim_index)
1660 {
1661   fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1662   fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1663   fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1664   fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1665   fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1666 }
1667 END_TEST
1668 */
1669 START_TEST(test_parse_list_null)
1670 {
1671   int num_param_args;
1672   char **param_args;
1673   parse_list_string(NULL, SEP, '{', '}',
1674                     &num_param_args, &param_args);
1675   fail_unless(num_param_args == 0, NULL);
1676   fail_unless(param_args == NULL, NULL);
1677 }
1678 END_TEST
1679 START_TEST(test_parse_list_single_simple)
1680 {
1681   int num_param_args;
1682   char **param_args;
1683   parse_list_string("arg", SEP, '{', '}',
1684                     &num_param_args, &param_args);
1685   fail_unless(num_param_args == 1, NULL);
1686   fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1687 }
1688 END_TEST
1689 START_TEST(test_parse_list_single_compound)
1690 {
1691   int num_param_args;
1692   char **param_args;
1693   parse_list_string("{x,y,z}", SEP, '{', '}',
1694                     &num_param_args, &param_args);
1695   fail_unless(num_param_args == 1, NULL);
1696   fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1697 }
1698 END_TEST
1699 START_TEST(test_parse_list_double_simple)
1700 {
1701   int num_param_args;
1702   char **param_args;
1703   parse_list_string("abc,def", SEP, '{', '}',
1704                     &num_param_args, &param_args);
1705   fail_unless(num_param_args == 2, NULL);
1706   fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1707   fail_unless(STRMATCH(param_args[1],"def"), NULL);
1708 }
1709 END_TEST
1710 @
1711 <<parse list string test case defs>>=
1712 TCase *tc_parse_list_string = tcase_create("parse list string");
1713 @
1714 <<parse list string test case add>>=
1715 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1716 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1717 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1718 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1719 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1720 suite_add_tcase(s, tc_parse_list_string);
1721 @
1722
1723
1724 \section{Unit tests}
1725
1726 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1727 <<checks>>=
1728 <<includes>>
1729 <<check includes>>
1730 <<definitions>>
1731 <<globals>>
1732 <<check globals>>
1733 <<check relative error>>
1734 <<functions>>
1735 <<test suite>>
1736 <<main check program>>
1737 @
1738
1739 <<check includes>>=
1740 #include <check.h>
1741 @
1742
1743 <<check globals>>=
1744 @
1745
1746 <<test suite>>=
1747 <<F tests>>
1748 <<determine dt tests>>
1749 <<happens tests>>
1750 <<does domain unfold tests>>
1751 <<randomly unfold domains tests>>
1752 <<suite definition>>
1753 @
1754
1755 <<suite definition>>=
1756 Suite *test_suite (void)
1757 {
1758   Suite *s = suite_create ("sawsim");
1759   <<F test case defs>>
1760   <<determine dt test case defs>>
1761   <<happens test case defs>>
1762   <<does domain unfold test case defs>>
1763   <<randomly unfold domains test case defs>>
1764
1765   <<F test case add>>
1766   <<determine dt test case add>>
1767   <<happens test case add>>
1768   <<does domain unfold test case add>>
1769   <<randomly unfold domains test case add>>
1770
1771 /*
1772   tcase_add_checked_fixture(tc_strip_address,
1773                             setup_strip_address,
1774                             teardown_strip_address);
1775 */
1776   return s;
1777 }
1778 @
1779
1780 <<main check program>>=
1781 int main(void)
1782 {
1783   int number_failed;
1784   Suite *s = test_suite();
1785   SRunner *sr = srunner_create(s);
1786   srunner_run_all(sr, CK_ENV);
1787   number_failed = srunner_ntests_failed(sr);
1788   srunner_free(sr);
1789   return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1790 }
1791 @
1792
1793 \subsection{F tests}
1794 <<F tests>>=
1795 <<worm-like chain tests>>
1796 @
1797 <<F test case defs>>=
1798 <<worm-like chain test case def>>
1799 @
1800 <<F test case add>>=
1801 <<worm-like chain test case add>>
1802 @
1803
1804 <<worm-like chain tests>>=
1805 extern double wlc(double x, double T, double p, double L);
1806 START_TEST(test_wlc_at_zero)
1807 {
1808   double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
1809   fail_unless(abs(wlc(x, T, p, L)) < lim, \
1810               "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
1811               x, T, p, L, abs(wlc(x, T, p, L)), lim);
1812 }
1813 END_TEST
1814
1815 START_TEST(test_wlc_at_half)
1816 {
1817   double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1818   /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1819    * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1820    *                = 0.25 * 3 = 3/4
1821    * linear term = x/L = 0.5
1822    * nonlinear + linear = 0.75 + 0.5 = 1.25
1823    * wlc = 10e21*1.25 = 12.5e21
1824    */
1825   fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1826               "wlc(%g, %g, %g, %g) = %g != %g",
1827                x, T, p, L, wlc(x, T, p, L), 12.5e21);
1828 }
1829 END_TEST
1830 @
1831 <<worm-like chain test case def>>=
1832 TCase *tc_wlc = tcase_create("WLC");
1833 @
1834
1835 <<worm-like chain test case add>>=
1836 tcase_add_test(tc_wlc, test_wlc_at_zero);
1837 tcase_add_test(tc_wlc, test_wlc_at_half);
1838 suite_add_tcase(s, tc_wlc);
1839 @
1840
1841 \subsection{Model tests}
1842
1843 Check the searching with [[linear_k]].
1844 Check overwhelming force treatment with the heavyside-esque step [[?]].
1845 <<determine dt tests>>=
1846 double linear_k(double F, environment_t *env, void *params)
1847 {
1848   double Fnot = *(double *)params;
1849   return Fnot+F;
1850 }
1851
1852 START_TEST(test_determine_dt_linear_k)
1853 {
1854   environment_t env;
1855   double dt_max=3.0, Fnot=3.0;
1856   double F[]={0,1,2,3,4,5,6};
1857   domain_t dom; /* use both parts at once for folded/unfolded */
1858   int i;
1859   env.T = 300.0;
1860 /*
1861   dom->next = dom->prev = NULL;
1862   dom->k_func_t = linear_k;
1863   dom->folded_params = &Fnot;
1864   dom->unfolded_params = !!!!!!!!!
1865   dom->destroy_folded = dom->destroy_unfolded = NULL;
1866   for( i=0; i < sizeof(F)/sizeof(double); i++) {
1867     fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1868   }
1869 */
1870 }
1871 END_TEST
1872 @
1873 <<determine dt test case defs>>=
1874 TCase *tc_determine_dt = tcase_create("Determine dt");
1875 @
1876 <<determine dt test case add>>=
1877 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1878 suite_add_tcase(s, tc_determine_dt);
1879 @
1880
1881 <<happens tests>>=
1882 @
1883 <<happens test case defs>>=
1884 @
1885 <<happens test case add>>=
1886 @
1887
1888 <<does domain unfold tests>>=
1889 @
1890 <<does domain unfold test case defs>>=
1891 @
1892 <<does domain unfold test case add>>=
1893 @
1894
1895 <<randomly unfold domains tests>>=
1896 @
1897 <<randomly unfold domains test case defs>>=
1898 @
1899 <<randomly unfold domains test case add>>=
1900 @
1901
1902
1903 \section{Balancing group extensions}
1904 \label{app.tension_balance}
1905
1906 For a given total extension $x$ (of the piezo), the various domain
1907 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
1908 amounts, and we need to tweak the portion that each extends to balance
1909 the tension amongst the active groups.  Since the tension balancing is
1910 separable from the bulk of the simulation, we break it out into a
1911 separate module.  The interface is defined in [[tension_balance.h]],
1912 the implementation in [[tension_balance.c]], and the unit testing in
1913 [[check_tension_balance.c]]
1914
1915 <<tension-balance.h>>=
1916 #ifndef TENSION_BALANCE_H
1917 #define TENSION_BALANCE_H
1918 <<license comment>>
1919 <<tension balance function declaration>>
1920 #endif /* TENSION_BALANCE_H */
1921 @
1922
1923 <<tension balance functions>>=
1924 <<one dimensional bracket>>
1925 <<one dimensional solve>>
1926 <<x of x0>>
1927 <<tension balance function>>
1928 @
1929
1930 <<tension balance module makefile lines>>=
1931 NW_SAWSIM_MODS += tension_balance
1932 CHECK_BINS += check_tension_balance
1933 check_tension_balance_MODS =
1934 @
1935
1936 The entire force balancing problem reduces to a solving two nested
1937 one-dimensional root-finding problems.  First we define the one
1938 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
1939 non-empty group).  $x(x_0)$ is also strictly monotonically increasing,
1940 so we can solve for $x_0$ such that $\sum_i x_i = x$.  [[last_x]]
1941 stores the last successful value of $x$, since we expect to be taking
1942 small steps ($x \approx$ [[last_x]]).  Setting [[last_x = -1]]
1943 indicates that the groups have changed.
1944 <<tension balance function declaration>>=
1945 double tension_balance(int num_tension_groups,
1946                        one_dim_fn_t **tension_handlers,
1947                        void **params,
1948                        double last_x,
1949                        double x);
1950 @
1951 <<tension balance function>>=
1952 double tension_balance(int num_tension_groups,
1953                        one_dim_fn_t **tension_handlers,
1954                        void **params,
1955                        double last_x,
1956                        double x)
1957 {
1958   static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
1959   double F, xo_guess, xo, lb, ub;
1960   double min_relx=1e-6, min_rely=1e-6;
1961   int max_steps=200, i;
1962
1963   assert(num_tension_groups > 0);
1964   assert(tension_handlers != NULL);
1965   assert(params != NULL);
1966   assert(num_tension_groups > 0);
1967
1968   if (num_tension_groups == 1) { /* only one group, no balancing required */
1969     xo = x;
1970   } else {
1971     //fprintf(stderr, "balancing for x = %g with ", x);
1972     //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1973     //fprintf(stderr, "\n");
1974     if (last_x == -1) { /* new group setup, reset x_xo_data */
1975       /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
1976       if (x_xo_data.xi != NULL)
1977         free(x_xo_data.xi);
1978       /* construct new x_xo_data */
1979       x_xo_data.n_groups = num_tension_groups;
1980       x_xo_data.tension_handler = tension_handlers;
1981       x_xo_data.handler_data = params;
1982       x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
1983       for (i=0; i<num_tension_groups; i++)
1984         x_xo_data.xi[i] = -1.0;
1985     }
1986     if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
1987       if (x == 0) xo_guess = length_scale;
1988       else        xo_guess = x/num_tension_groups;
1989 #ifdef DEBUG
1990       fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
1991 #endif
1992       oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
1993     } else { /* work off the last balanced point */
1994 #ifdef DEBUG
1995       fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
1996 #endif
1997       if (x > last_x) {
1998         lb = x_xo_data.xi[0];
1999         ub = x_xo_data.xi[0]+ x-last_x;   /* apply all change to x0 */
2000       } else if (x < last_x) {
2001         lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2002         ub = x_xo_data.xi[0];
2003       } else { /* x == last_x */
2004         fprintf(stderr, "not moving\n");
2005         lb= x_xo_data.xi[0];
2006         ub= x_xo_data.xi[0];
2007       }
2008     }
2009 #ifdef DEBUG
2010     fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2011             __FUNCTION__, x, lb, ub);
2012 #endif
2013     xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2014     F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2015 #ifdef DEBUG
2016     if (num_tension_groups > 1) {
2017       fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2018       for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2019       fprintf(stderr, "\n");
2020     }
2021 #endif
2022   }
2023   F = (*tension_handlers[0])(xo, params[0]);
2024   return F;
2025 }
2026 @
2027
2028 <<tension balance internal definitions>>=
2029 <<x of x0 definitions>>
2030 @
2031
2032 <<x of x0 definitions>>=
2033 typedef struct x_of_xo_data_struct {
2034   int n_groups;
2035   one_dim_fn_t **tension_handler; /* array of fn pointers */
2036   void **handler_data;            /* array of void* pointers */
2037   double *xi;
2038 } x_of_xo_data_t;
2039 @
2040
2041 <<x of x0>>=
2042 double x_of_xo(double xo, void *pdata)
2043 {
2044   x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2045   double F, x=0, xi, xi_guess, lb, ub;
2046   int i;
2047   double min_relx=1e-6, min_rely=1e-6;
2048   int max_steps=200;
2049   assert(data->n_groups > 0);
2050   data->xi[0] = xo;
2051   F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2052 #ifdef DEBUG
2053   fprintf(stderr, "%s: finding x(x0=%g).  F0(x0) = %g\n", __FUNCTION__, xo, F);
2054 #endif
2055   x += xo;
2056   if (data->xi)  data->xi[0] = xo;
2057   for (i=1; i < data->n_groups; i++) {
2058     if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2059     else                   xi_guess = data->xi[i];
2060 #ifdef DEBUG
2061   fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
2062 #endif
2063     oneD_bracket(data->tension_handler[i], data->handler_data[i], F,  xi_guess, &lb, &ub);
2064     xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2065 #ifdef DEBUG
2066   fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2067 #endif
2068     data->xi[i] = xi;
2069     x += xi;
2070     if (data->xi)  data->xi[i] = xi;
2071   }
2072 #ifdef DEBUG
2073   fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2074 #endif
2075   return x;
2076 }
2077 @
2078
2079 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2080 number of steps for monotonically increasing $f(x)$.  Simple
2081 bisection, so it's robust and converges fairly quickly.  You can set
2082 the maximum number of steps to take, but if you have enough processor
2083 power, it's better to limit your precision with the relative limits
2084 [[min_relx]] and [[y]].  These ensure that the width of the bracket is
2085 small on the length scale of it's larger side.  Note that \emph{both}
2086 relative limits must be satisfied for the search to stop.
2087 <<one dimensional function definition>>=
2088 /* equivalent to gsl_func_t */
2089 typedef double one_dim_fn_t(double x, void *params);
2090 @
2091 <<one dimensional solve declaration>>=
2092 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2093                   double min_relx, double min_rely, int max_steps,
2094                   double *pYx);
2095 @
2096
2097 <<one dimensional solve>>=
2098 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2099 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2100                   double min_relx, double min_rely, int max_steps,
2101                   double *pYx)
2102 {
2103   double yx, ylb, yub, x;
2104   int i=0;
2105   assert(ub >= lb);
2106   ylb = (*f)(lb, params);
2107   yub = (*f)(ub, params);
2108
2109   /* check some simple cases */
2110   if (ylb == yub) {
2111     if (ylb != y)  return HUGE_VAL; /* error! f(x)=y not bounded */
2112     else           return lb; /* any x on the range [lb, ub] would work */
2113   }
2114   if (ylb == y) { x=lb; yx=ylb; goto end; }
2115   if (yub == y) { x=ub; yx=yub; goto end; }
2116
2117 #ifdef DEBUG
2118   fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2119 #endif
2120   assert(ylb < y); assert(yub > y);
2121   for (i=0; i<max_steps; i++) {
2122     x = (lb+ub)/2.0;
2123     yx = (*f)(x, params);
2124 #ifdef DEBUG
2125   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);
2126 #endif
2127     if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2128     else if (yx > y)  { ub=x; yub=yx; }
2129     else /* yx < y */ { lb=x; ylb=yx; }
2130   }
2131  end:
2132   if (pYx) *pYx = yx;
2133 #ifdef DEBUG
2134   fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2135 #endif
2136   return x;
2137 }
2138 @
2139
2140 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2141 Generate bracketing $x$ values through bisection or exponential growth.
2142 <<one dimensional bracket declaration>>=
2143 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2144 @
2145
2146 <<one dimensional bracket>>=
2147 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2148 {
2149   double yx, step, x=xguess;
2150   int i=0;
2151   yx = (*f)(x, params);
2152 #ifdef DEBUG
2153   fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2154 #endif
2155   if (yx > y) {
2156     assert(x > 0.0);
2157     *ub = x;
2158     *lb = 0;
2159 #ifdef DEBUG
2160     fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2161 #endif
2162   } else {
2163     *lb = x;
2164     if (x == 0) x = length_scale/2.0; /* HACK */
2165     while (yx < y) {
2166       x *= 2.0;
2167       yx = (*f)(x, params);
2168 #ifdef DEBUG
2169       fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2170 #endif
2171     }
2172     *ub = x;
2173 #ifdef DEBUG
2174     fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2175 #endif
2176   }
2177   //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2178 }
2179 @
2180
2181 \subsection{Balancing implementation}
2182
2183 <<tension-balance.c>>=
2184 //#define DEBUG
2185 <<license comment>>
2186 <<tension balance includes>>
2187 #include "tension_balance.h"
2188
2189 double length_scale = 1e-10; /* HACK */
2190
2191 <<tension balance internal definitions>>
2192 <<tension balance functions>>
2193 @
2194
2195 <<tension balance includes>>=
2196 #include <assert.h> /* assert()          */
2197 #include <stdlib.h> /* NULL              */
2198 #include <math.h>   /* HUGE_VAL, macro constant, so don't need to link to math */
2199 #include <stdio.h>  /* fprintf(), stdout */
2200 #include "global.h"
2201 @
2202
2203 \subsection{Balancing unit tests}
2204
2205 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2206 <<check-tension-balance.c>>=
2207 <<tension balance check includes>>
2208 <<tension balance test suite>>
2209 <<main check program>>
2210 @
2211
2212 <<tension balance check includes>>=
2213 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2214 #include <assert.h> /* assert()                      */
2215 <<check includes>>
2216 #include "global.h"
2217 #include "tension_balance.h"
2218 @
2219
2220 <<tension balance test suite>>=
2221 <<tension balance function tests>>
2222 <<tension balance suite definition>>
2223 @
2224
2225 <<tension balance suite definition>>=
2226 Suite *test_suite (void)
2227 {
2228   Suite *s = suite_create ("tension balance");
2229   <<tension balance function test case defs>>
2230
2231   <<tension balance function test case adds>>
2232   return s;
2233 }
2234 @
2235
2236 <<tension balance function tests>>=
2237 <<check relative error>>
2238
2239 double hooke(double x, void *pK)
2240 {
2241   assert(x >= 0);
2242   return *((double*)pK) * x;
2243 }
2244
2245 START_TEST(test_single_function)
2246 {
2247   double k=5, x=3, last_x=2.0, F;
2248   one_dim_fn_t *handlers[] = {&hooke};
2249   void *data[] =             {&k};
2250   F = tension_balance(1, handlers, data, last_x, x);
2251   fail_unless(F = k*x, NULL);
2252 }
2253 END_TEST
2254 @
2255
2256 We can also test balancing two springs with different spring constants.
2257 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2258 <<tension balance function tests>>=
2259 START_TEST(test_double_hooke)
2260 {
2261   double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2262   one_dim_fn_t *handlers[] = {&hooke, &hooke};
2263   void *data[] =             {&k1,    &k2};
2264   F = tension_balance(2, handlers, data, last_x, x);
2265   x1e = x*k2/(k1+k2);
2266   Fe = k1*x1e;
2267   x2e = x1e*k1/k2;
2268   //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2269   //CHECK_ERR(1e-6, x1e, xi[0]);
2270   //CHECK_ERR(1e-6, x2e, xi[1]);
2271   CHECK_ERR(1e-6, Fe, F);
2272 }
2273 END_TEST
2274 @
2275 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2276
2277 <<tension balance function test case defs>>=
2278 TCase *tc_tbfunc = tcase_create("tension balance function");
2279 @
2280
2281 <<tension balance function test case adds>>=
2282 tcase_add_test(tc_tbfunc, test_single_function);
2283 tcase_add_test(tc_tbfunc, test_double_hooke);
2284 suite_add_tcase(s, tc_tbfunc);
2285 @
2286
2287 \section{Lists}
2288
2289 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2290 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2291 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2292
2293 <<list.h>>=
2294 #ifndef LIST_H
2295 #define LIST_H
2296 <<license comment>>
2297 <<list definitions>>
2298 <<list declarations>>
2299 #endif /* LIST_H */
2300 @
2301
2302 <<list declarations>>=
2303 <<head and tail declarations>>
2304 <<list length declaration>>
2305 <<push declaration>>
2306 <<pop declaration>>
2307 <<list destroy declaration>>
2308 <<list to array declaration>>
2309 <<move declaration>>
2310 <<sort declaration>>
2311 <<unique declaration>>
2312 <<list print declaration>>
2313 @
2314
2315 <<list functions>>=
2316 <<create and destroy node>>
2317 <<head and tail>>
2318 <<list length>>
2319 <<push>>
2320 <<pop>>
2321 <<list destroy>>
2322 <<list to array>>
2323 <<move>>
2324 <<sort>>
2325 <<unique>>
2326 <<list print>>
2327 @
2328
2329 <<list module makefile lines>>=
2330 NW_SAWSIM_MODS += list
2331 CHECK_BINS += check_list
2332 check_list_MODS =
2333 @
2334
2335 <<list definitions>>=
2336 typedef struct list_struct {
2337   struct list_struct *next;
2338   struct list_struct *prev;
2339   void *d; /* data */
2340 } list_t;
2341
2342 <<comparison function typedef>>
2343 @
2344
2345 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2346 <<head and tail declarations>>=
2347 list_t *head(list_t *list);
2348 list_t *tail(list_t *list);
2349 @
2350 <<head and tail>>=
2351 list_t *head(list_t *list)
2352 {
2353   if (list == NULL)
2354     return list;
2355   while (list->prev != NULL) {
2356     list = list->prev;
2357   }
2358   return list;
2359 }
2360
2361 list_t *tail(list_t *list)
2362 {
2363   if (list == NULL)
2364     return list;
2365   while (list->next != NULL) {
2366     list = list->next;
2367   }
2368   return list;
2369 }
2370 @
2371
2372 <<list length declaration>>=
2373 int list_length(list_t *list);
2374 @
2375 <<list length>>=
2376 int list_length(list_t *list)
2377 {
2378   int i;
2379   if (list == NULL)
2380     return 0;
2381   list = head(list);
2382   i = 1;
2383   while (list->next != NULL) {
2384     list = list->next;
2385     i += 1;
2386   }
2387   return i;
2388 }
2389 @
2390
2391 [[push]] inserts a node relative to the current position in the list
2392 without changing the current position.
2393 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2394 so we set it to that of the pushed domain.
2395 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2396 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2397 <<push declaration>>=
2398 void push(list_t **pList, void *data, int below);
2399 @
2400 <<push>>=
2401 void push(list_t **pList, void *data, int below)
2402 {
2403   list_t *list, *new_node;
2404   assert(pList != NULL);
2405   list = *pList;
2406   new_node = create_node(data);
2407   if (list == NULL)
2408     *pList = new_node;
2409   else if (below==0) { /* insert above */
2410     if (list->prev != NULL)
2411       list->prev->next = new_node;
2412     new_node->prev = list->prev;
2413     list->prev = new_node;
2414     new_node->next = list;
2415   } else {         /* insert below */
2416     if (list->next != NULL)
2417       list->next->prev = new_node;
2418     new_node->next = list->next;
2419     list->next = new_node;
2420     new_node->prev = list;
2421   }
2422 }
2423 @
2424
2425 [[pop]] removes the current domain node, moving the current position
2426 to the node after it, unless that node is [[NULL]], in which case move
2427 the current position to the node before it.
2428 <<pop declaration>>=
2429 void *pop(list_t **pList);
2430 @
2431 <<pop>>=
2432 void *pop(list_t **pList)
2433 {
2434   list_t *list, *popped;
2435   void *data;
2436   assert(pList != NULL);
2437   list = *pList;
2438   assert(list != NULL); /* not an empty list */
2439   popped = list;
2440   /* bypass the popped node */
2441   if (list->prev != NULL)
2442      list->prev->next = list->next;
2443   if (list->next != NULL) {
2444      list->next->prev = list->prev;
2445      *pList = list->next;
2446   } else
2447      *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2448   data = popped->d;
2449   destroy_node(popped);
2450   return data;
2451 }
2452 @
2453
2454 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2455 If the cleanup function is [[NULL]], no cleanup function is called.
2456 The nodes are not popped in any particular order.
2457 <<list destroy declaration>>=
2458 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2459 @
2460 <<list destroy>>=
2461 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2462 {
2463   list_t *list;
2464   void *data;
2465   assert(pList != NULL);
2466   list = *pList;
2467   *pList = NULL;
2468   assert(list != NULL); /* not an empty list */
2469   while (list != NULL) {
2470     data = pop(&list);
2471     if (destroy != NULL)
2472       destroy(data);
2473   }
2474 }
2475 @
2476
2477 [[list_to_array]] converts a list to an array of pointers, preserving order.
2478 <<list to array declaration>>=
2479 void list_to_array(list_t *list, int *length, void ***array);
2480 @
2481 <<list to array>>=
2482 void list_to_array(list_t *list, int *pLength, void ***pArray)
2483 {
2484   void **array;
2485   int i,length;
2486   assert(list != NULL);
2487   assert(pLength != NULL);
2488   assert(pArray != NULL);
2489
2490   length = list_length(list);
2491   array = (void **)malloc(sizeof(void *)*length);
2492   assert(array != NULL);
2493   list = head(list);
2494   i=0;
2495   while (list != NULL) {
2496     array[i++] = list->d;
2497     list = list->next;
2498   }
2499   *pLength = length;
2500   *pArray = array;
2501 }
2502 @
2503
2504 [[move]] moves the current element along the list in either direction.
2505 <<move declaration>>=
2506 void move(list_t **pList, int downstream);
2507 @
2508 <<move>>=
2509 void move(list_t **pList, int downstream)
2510 {
2511   list_t *list, *flipper;
2512   void *data;
2513   assert(pList != NULL);
2514   list = *pList;          /* a>B>c>d */
2515   assert(list != NULL); /* not an empty list */
2516   if (downstream == 0)
2517     flipper = list->prev; /* flipper is A */
2518   else
2519     flipper = list->next; /* flipper is C */
2520   /* ensure there is somewhere to go in the direction we're heading */
2521   assert(flipper != NULL);
2522   /* remove the current element */
2523   data = pop(&list);      /* data=B, a>C>d */
2524   /* and put it back in in the correct spot */
2525   list = flipper;
2526   if (downstream == 0) {
2527     push(&list, data, 0); /* b>A>c>d */
2528     list = list->prev;    /* B>a>c>d   */
2529   } else {
2530     push(&list, data, 1); /* a>C>b>d */
2531     list = list->next;    /* a>c>B>d */
2532   }
2533   *pList = list;
2534 }
2535 @
2536
2537 [[sort]] sorts a list from smallest to largest according to a user
2538 specified comparison.
2539 <<comparison function typedef>>=
2540 typedef int (comparison_fn_t)(void *A, void *B);
2541 @
2542 For example
2543 <<int comparison function>>=
2544 int int_comparison(void *A, void *B)
2545 {
2546   int a,b;
2547   a = *((int *)A);
2548   b = *((int *)B);
2549   if (a > b) return 1;
2550   else if (a == b) return 0;
2551   else return -1;
2552 }
2553 @
2554 <<sort declaration>>=
2555 void sort(list_t **pList, comparison_fn_t *comp);
2556 @
2557 Here we use the bubble sort algorithm.
2558 <<sort>>=
2559 void sort(list_t **pList, comparison_fn_t *comp)
2560 {
2561   list_t *list;
2562   int stable=0;
2563   assert(pList != NULL);
2564   list = *pList;
2565   assert(list != NULL); /* not an empty list */
2566   while (stable == 0) {
2567     stable = 1;
2568     list = head(list);
2569     while (list->next != NULL) {
2570       if ((*comp)(list->d, list->next->d) > 0) {
2571         stable = 0;
2572         move(&list, 1 /* downstream */);
2573       } else
2574         list = list->next;
2575     }
2576   }
2577   *pList = list;
2578 }
2579 @
2580
2581 [[unique]] removes duplicates from a sorted list according to a user
2582 specified comparison.  Don't do this unless you have other pointers
2583 any dynamically allocated data in your list, because [[unique]] just
2584 drops any non-unique data without freeing it.
2585 <<unique declaration>>=
2586 void unique(list_t **pList, comparison_fn_t *comp);
2587 @
2588 <<unique>>=
2589 void unique(list_t **pList, comparison_fn_t *comp)
2590 {
2591   list_t *list;
2592   assert(pList != NULL);
2593   list = *pList;
2594   assert(list != NULL); /* not an empty list */
2595   list = head(list);
2596   while (list->next != NULL) {
2597     if ((*comp)(list->d, list->next->d) == 0) {
2598       pop(&list);
2599     } else
2600       list = list->next;
2601   }
2602   *pList = list;
2603 }
2604 @
2605
2606 [[list_print]] prints a list to a [[FILE *]] stream.
2607 <<list print declaration>>=
2608 void list_print(FILE *file, list_t *list, char *list_name);
2609 @
2610 <<list print>>=
2611 void list_print(FILE *file, list_t *list, char *list_name)
2612 {
2613   fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2614   list = head(list);
2615   while (list != NULL) {
2616     fprintf(file, " %p", list->d);
2617     list = list->next;
2618   }
2619   fprintf(file, "\n");
2620 }
2621 @
2622
2623 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2624 <<create and destroy node>>=
2625 list_t *create_node(void *data)
2626 {
2627   list_t *ret = (list_t *)malloc(sizeof(list_t));
2628   assert(ret != NULL);
2629   ret->prev = NULL;
2630   ret->next = NULL;
2631   ret->d = data;
2632   return ret;
2633 }
2634
2635 void destroy_node(list_t *node)
2636 {
2637   if (node)
2638     free(node);
2639 }
2640 @
2641 The user must free the data pointed to by the node on their own.
2642
2643 \subsection{List implementation}
2644
2645 <<list.c>>=
2646 <<license comment>>
2647 <<list includes>>
2648 #include "list.h"
2649 <<list functions>>
2650 @
2651
2652 <<list includes>>=
2653 #include <assert.h> /* assert()                                */
2654 #include <stdlib.h> /* malloc(), free(), rand()                */
2655 #include <stdio.h>  /* fprintf(), stdout, FILE                 */
2656 #include "global.h" /* destroy_data_func_t */
2657 @
2658
2659 \subsection{List unit tests}
2660
2661 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2662 <<check-list.c>>=
2663 <<list check includes>>
2664 <<list test suite>>
2665 <<main check program>>
2666 @
2667
2668 <<list check includes>>=
2669 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2670 #include <stdio.h>  /* FILE                          */
2671 <<check includes>>
2672 #include "global.h"
2673 #include "list.h"
2674 @
2675
2676 <<list test suite>>=
2677 <<push tests>>
2678 <<pop tests>>
2679 <<list suite definition>>
2680 @
2681
2682 <<list suite definition>>=
2683 Suite *test_suite (void)
2684 {
2685   Suite *s = suite_create ("list");
2686   <<push test case defs>>
2687   <<pop test case defs>>
2688
2689   <<push test case adds>>
2690   <<pop test case adds>>
2691   return s;
2692 }
2693 @
2694
2695 <<push tests>>=
2696 START_TEST(test_push)
2697 {
2698   list_t *list=NULL;
2699   int i, p, e, n=3;
2700   for (i=0; i<n; i++)
2701     push(&list, (void *)i, 1);
2702   fail_unless(list == head(list), NULL);
2703   fail_unless( (int)list->d == 0 );
2704   for (i=0; i<n; i++) {
2705     /* we expect 0, n-1, n-2, ..., 1 */
2706     if (i == 0) e = 0;
2707     else        e = n-i;
2708     fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2709   }
2710 }
2711 END_TEST
2712 @
2713
2714 <<push test case defs>>=
2715 TCase *tc_push = tcase_create("push");
2716 @
2717
2718 <<push test case adds>>=
2719 tcase_add_test(tc_push, test_push);
2720 suite_add_tcase(s, tc_push);
2721 @
2722
2723 <<pop tests>>=
2724 @
2725 <<pop test case defs>>=
2726 @
2727 <<pop test case adds>>=
2728 @
2729
2730 \section{Function string evaluation}
2731
2732 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).
2733 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2734 The solution is to run a scripting language as a subprocess accessed via pipes.
2735 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2736
2737 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2738 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2739 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.
2740 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2741
2742 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.
2743 Then you can either statically or dynamically link to those hardcoded functions.
2744 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2745
2746 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2747 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2748
2749 <<string-eval.h>>=
2750 #ifndef STRING_EVAL_H
2751 #define STRING_EVAL_H
2752 <<license comment>>
2753 <<string eval setup declaration>>
2754 <<string eval function declaration>>
2755 <<string eval teardown declaration>>
2756 #endif /* STRING_EVAL_H */
2757 @
2758
2759 <<string eval module makefile lines>>=
2760 NW_SAWSIM_MODS += string_eval
2761 CHECK_BINS += check_string_eval
2762 check_string_eval_MODS =
2763 @
2764
2765 For an introduction to POSIX process control, see\\
2766  \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2767  \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2768  and of course, the relavent [[man]] pages.
2769
2770 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2771 [[execvp]] replaces the calling process' program with a new program.
2772 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2773 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2774  but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2775 The new program needs command line arguments, just like it would if you were running it from a shell.
2776 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2777 with the final array entry being a [[NULL]] pointer.
2778
2779 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2780 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2781 <<string eval subprocess definitions>>=
2782 #define SUBPROCESS "python"
2783 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2784 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2785 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2786 @
2787 The [[i]] option lets Python know that it should run in interactive mode.
2788 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2789 In interactive mode, python acts on every instruction as soon as it is recieved.
2790 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.
2791 %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.
2792
2793 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2794 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2795 [[fork]] returns two copies of the same program, executing the original code.
2796 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2797 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2798
2799 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.
2800 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2801 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2802 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2803 <<string eval pipe definitions>>=
2804 #define PIPE_READ  0   /* the end you read from */
2805 #define PIPE_WRITE 1   /* the end you write to */
2806
2807 #define STDIN      0   /* initial index of stdin pair  */
2808 #define STDOUT     2   /* initial index of stdout pair */
2809
2810 @
2811 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2812
2813 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.
2814
2815 <<string eval setup declaration>>=
2816 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2817 @
2818 <<string eval setup definition>>=
2819 void string_eval_setup(FILE **pIn, FILE **pOut)
2820 {
2821   pid_t pid;
2822   int pfd[4]; /* file descriptors for stdin and stdout */
2823   int rc;
2824   assert(pipe(pfd+STDIN)  != -1); /* stdin pair  (in, out) */
2825   assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2826
2827   pid = fork(); /* split process into two copies */
2828
2829   if (pid == -1) { /* fork error */
2830     perror("fork error");
2831     exit(1);
2832   } else if (pid == 0) { /* child */
2833     close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input   */
2834     close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2835     dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin  (closes original stdin) */
2836     dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2837     execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2838     perror("exec error"); /* exec shouldn't return */
2839     _exit(1);
2840   } else { /* parent */
2841     close(pfd[STDIN+PIPE_READ]);   /* close stdin pipe output */
2842     close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2843     *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2844     if ( *pIn == NULL ) {
2845       perror("fdopen (in)");
2846       exit(1);
2847     }
2848     *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2849     if ( *pOut == NULL ) {
2850       perror("fdopen (out)");
2851       exit(1);
2852     }
2853   }
2854 }
2855 @
2856
2857 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2858 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2859 <<string eval function declaration>>=
2860 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2861 @
2862 <<string eval function definition>>=
2863 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2864 {
2865   int rc;
2866   rc = fprintf(in, "%s", input);
2867   assert(rc == strlen(input));
2868   fflush(in);
2869   fflush(out);
2870   alarm(1); /* set a one second timeout on the read */
2871   assert( fgets(output, buflen, out) != NULL );
2872   alarm(0); /* cancel the timeout */
2873   //fprintf(stderr, "eval: %s ----> %s", input, output);
2874 }
2875 @
2876 The [[alarm]] calls set and clear a timeout on the returned output.
2877 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.
2878 This protects against invalid input for which a line of output is not printed to [[stdout]].
2879 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2880 If you are getting strange results, check your python code seperately. TODO, better error handling.
2881
2882 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2883 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2884 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2885 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2886 <<string eval teardown declaration>>=
2887 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2888 @
2889
2890 <<string eval teardown definition>>=
2891 void string_eval_teardown(FILE **pIn, FILE **pOut)
2892 {
2893   pid_t pid=0;
2894   int stat_loc;
2895
2896   /* redirect Python's stderr */
2897   fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2898   fflush(*pIn);
2899
2900   /* close pipes */
2901   assert( fclose(*pIn) == 0 );
2902   *pIn = NULL;
2903   assert( fclose(*pOut) == 0 );
2904   *pOut = NULL;
2905
2906   /* wait for python to exit */
2907   while (pid <= 0) {
2908     pid = wait(&stat_loc);
2909     if (pid < 0) {
2910       perror("pid");
2911     }
2912   }
2913
2914   /*
2915   if (WIFEXITED(stat_loc)) {
2916     printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2917   } else if (WIFSIGNALED(stat_loc)) {
2918     printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2919   }
2920   */
2921 }
2922 @
2923 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2924
2925 \subsection{String evaluation implementation}
2926
2927 <<string-eval.c>>=
2928 <<license comment>>
2929 <<string eval includes>>
2930 #include "string_eval.h"
2931 <<string eval internal definitions>>
2932 <<string eval functions>>
2933 @
2934
2935 <<string eval includes>>=
2936 #include <assert.h> /* assert()                    */
2937 #include <stdlib.h> /* NULL                        */
2938 #include <stdio.h>  /* fprintf(), stdout, fdopen() */
2939 #include <string.h> /* strlen()                    */
2940 #include <sys/types.h> /* pid_t                    */
2941 #include <unistd.h> /* pipe(), fork(), execvp(), alarm()    */
2942 #include <wait.h>   /* wait()                      */
2943 @
2944
2945 <<string eval internal definitions>>=
2946 <<string eval subprocess definitions>>
2947 <<string eval pipe definitions>>
2948 @
2949
2950 <<string eval functions>>=
2951 <<string eval setup definition>>
2952 <<string eval function definition>>
2953 <<string eval teardown definition>>
2954 @
2955
2956 \subsection{String evaluation unit tests}
2957
2958 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2959 <<check-string-eval.c>>=
2960 <<string eval check includes>>
2961 <<string comparison definition and globals>>
2962 <<string eval test suite>>
2963 <<main check program>>
2964 @
2965
2966 <<string eval check includes>>=
2967 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2968 #include <stdio.h>  /* printf()                              */
2969 #include <assert.h> /* assert()                              */
2970 #include <string.h> /* strlen()                              */
2971 #include <signal.h> /* SIGKILL                               */
2972 <<check includes>>
2973 #include "string_eval.h"
2974 @
2975
2976 <<string eval test suite>>=
2977 <<string eval tests>>
2978 <<string eval suite definition>>
2979 @
2980
2981 <<string eval suite definition>>=
2982 Suite *test_suite (void)
2983 {
2984   Suite *s = suite_create ("string eval");
2985   <<string eval test case defs>>
2986
2987   <<string eval test case add>>
2988   return s;
2989 }
2990 @
2991
2992 <<string eval tests>>=
2993 START_TEST(test_setup_teardown)
2994 {
2995   FILE *in, *out;
2996   string_eval_setup(&in, &out);
2997   string_eval_teardown(&in, &out);
2998 }
2999 END_TEST
3000 START_TEST(test_invalid_command)
3001 {
3002   FILE *in, *out;
3003   char input[80], output[80]={};
3004   string_eval_setup(&in, &out);
3005   sprintf(input, "print ABCDefg\n");
3006   string_eval(in, out, input, 80, output);
3007   string_eval_teardown(&in, &out);
3008 }
3009 END_TEST
3010 START_TEST(test_simple_eval)
3011 {
3012   FILE *in, *out;
3013   char input[80], output[80]={};
3014   string_eval_setup(&in, &out);
3015   sprintf(input, "print 3+4*5\n");
3016   string_eval(in, out, input, 80, output);
3017   fail_unless(STRMATCH(output,"23\n"), NULL);
3018   string_eval_teardown(&in, &out);
3019 }
3020 END_TEST
3021 START_TEST(test_multiple_evals)
3022 {
3023   FILE *in, *out;
3024   char input[80], output[80]={};
3025   string_eval_setup(&in, &out);
3026   sprintf(input, "print 3+4*5\n");
3027   string_eval(in, out, input, 80, output);
3028   fail_unless(STRMATCH(output,"23\n"), NULL);
3029   sprintf(input, "print (3**2 + 4**2)**0.5\n");
3030   string_eval(in, out, input, 80, output);
3031   fail_unless(STRMATCH(output,"5.0\n"), NULL);
3032   string_eval_teardown(&in, &out);
3033 }
3034 END_TEST
3035 START_TEST(test_eval_with_state)
3036 {
3037   FILE *in, *out;
3038   char input[80], output[80]={};
3039   string_eval_setup(&in, &out);
3040   sprintf(input, "print 3+4*5\n");
3041   fprintf(in, "A = 3\n");
3042   sprintf(input, "print A*3\n");
3043   string_eval(in, out, input, 80, output);
3044   fail_unless(STRMATCH(output,"9\n"), NULL);
3045   string_eval_teardown(&in, &out);
3046 }
3047 END_TEST
3048 @
3049 <<string eval test case defs>>=
3050 TCase *tc_string_eval = tcase_create("string_eval");
3051 @
3052 <<string eval test case add>>=
3053 tcase_add_test(tc_string_eval, test_setup_teardown);
3054 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3055 tcase_add_test(tc_string_eval, test_simple_eval);
3056 tcase_add_test(tc_string_eval, test_multiple_evals);
3057 tcase_add_test(tc_string_eval, test_eval_with_state);
3058 suite_add_tcase(s, tc_string_eval);
3059 @
3060
3061
3062 \section{Accelerating function evaluation}
3063
3064 My first version-0.3 code was running very slowly.
3065 With the options suggested in the help
3066 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3067 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3068 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3069 That's only 15 calls per solution, so the search algorithm seems reasonable.
3070 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3071
3072 <<accel-k.h>>=
3073 #ifndef ACCEL_K_H
3074 #define ACCEL_K_H
3075 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3076 void free_accels();
3077 #endif /* ACCEL_K_H */
3078 @
3079
3080 <<accel k module makefile lines>>=
3081 NW_SAWSIM_MODS += accel_k
3082 #CHECK_BINS += check_accel_k
3083 check_accel_k_MODS =
3084 @
3085
3086 <<accel-k.c>>=
3087 #include <assert.h>  /* assert()                */
3088 #include <stdlib.h>  /* realloc(), free(), NULL */
3089 #include "global.h"  /* environment_t           */
3090 #include "k_model.h" /* k_func_t                */
3091 #include "interp.h"  /* interp_*                */
3092 #include "accel_k.h"
3093
3094 typedef struct accel_k_struct {
3095   interp_table_t *itable;
3096   k_func_t *k;
3097   environment_t *env;
3098   void *params;
3099 } accel_k_t;
3100
3101 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3102 static int num_accels = 0;
3103 static accel_k_t *accels=NULL;
3104
3105 /* Wrap k in the standard f(x) acceleration form */
3106 static double k_wrap(double F, void *params)
3107 {
3108   accel_k_t *p = (accel_k_t *)params;
3109   return (*p->k)(F, p->env, p->params);
3110 }
3111
3112 static int k_tol(double FA, double kA, double FB, double kB)
3113 {
3114   assert(FB > FA);
3115   if (FB-FA > 1e-12) {
3116     //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3117     return 1; /* unnacceptable */
3118   } else {
3119     //printf("acceptable tol\n");
3120     return 0; /* acceptable */
3121   }
3122 }
3123
3124 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3125 {
3126   int i=num_accels;
3127   accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3128   assert(accels != NULL);
3129   accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3130   accels[i].k = k;
3131   accels[i].env = env;
3132   accels[i].params = params;
3133   return i;
3134 }
3135
3136 void free_accels()
3137 {
3138   int i;
3139   for (i=0; i<num_accels; i++)
3140     interp_table_free(accels[i].itable);
3141   num_accels=0;
3142   if (accels) free(accels);
3143   accels = NULL;
3144 }
3145
3146 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3147 {
3148   int i;
3149   for (i=0; i<num_accels; i++) {
3150     if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3151       /* We've already setup this function.
3152        * WARNING: we're assuming param and env are the same. */
3153       return interp_table_eval(accels[i].itable, accels+i, F);
3154     }
3155   }
3156   /* set up a new acceleration environment */
3157   i = add_accel_k(k, env, params);
3158   return interp_table_eval(accels[i].itable, accels+i, F);
3159 }
3160 @
3161
3162 \section{Tension models}
3163 \label{sec.tension_models}
3164
3165 TODO: tension model intro.
3166 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.
3167
3168 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3169 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]].
3170
3171 <<tension-model.h>>=
3172 #ifndef TENSION_MODEL_H
3173 #define TENSION_MODEL_H
3174 <<license comment>>
3175 <<tension handler types>>
3176 <<constant tension model declarations>>
3177 <<hooke tension model declarations>>
3178 <<worm-like chain tension model declarations>>
3179 <<freely-jointed chain tension model declarations>>
3180 <<find tension definitions>>
3181 <<tension model global declarations>>
3182 #endif /* TENSION_MODEL_H */
3183 @
3184
3185 <<tension model module makefile lines>>=
3186 NW_SAWSIM_MODS += tension_model
3187 #CHECK_BINS += check_tension_model
3188 check_tension_model_MODS =
3189 @
3190 <<tension model utils makefile lines>>=
3191 TENSION_MODEL_MODS = tension_model parse list tension_balance
3192 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3193         $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3194         $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3195 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3196         notangle -Rtension-model-utils.c $< > $@
3197 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3198         gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3199 clean_tension_model_utils :
3200         rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3201 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3202 case, the directories) as ‘order-only’ prerequisites.  The timestamp
3203 on these prerequisits does not effect whether the rules are executed.
3204 This is appropriate for directories, where we don't need to recompile
3205 after an unrelated has been added to the directory, but only when the
3206 source prerequisites change.  See the [[make]] documentation for more
3207 details
3208 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3209
3210 \subsection{Null}
3211 \label{sec.null_tension}
3212
3213 For unstretchable domains.
3214
3215 <<null tension model getopt>>=
3216 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3217 @
3218
3219 \subsection{Constant}
3220 \label{sec.const_tension}
3221
3222 <<constant tension functions>>=
3223 <<constant tension function>>
3224 <<constant tension parameter create/destroy functions>>
3225 @
3226
3227 <<constant tension model declarations>>=
3228 <<constant tension function declaration>>
3229 <<constant tension parameter create/destroy function declarations>>
3230 <<constant tension model global declarations>>
3231
3232 @
3233
3234 An infinitely stretchable domain providing a constant tension.
3235
3236 <<constant tension function declaration>>=
3237 extern double const_tension_handler(double x, void *pdata);
3238 @
3239 <<constant tension function>>=
3240 double const_tension_handler(double x, void *pdata)
3241 {
3242   tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3243   list_t *list = data->group;
3244   double F;
3245
3246   assert (x >= 0.0);
3247   list = head(list);
3248   assert(list != NULL); /* empty active group?! */
3249   F = ((const_tension_param_t *)list->d)->F;
3250 #ifdef DEBUG
3251   fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3252 #endif
3253   while (list != NULL) {
3254     assert(((const_tension_param_t *)list->d)->F == F);
3255     list = list->next;
3256   }
3257   return F;
3258 }
3259
3260 @
3261
3262 <<constant tension structure>>=
3263 typedef struct const_tension_param_struct {
3264   double F; /* tension (force) in N */
3265 } const_tension_param_t;
3266 @
3267
3268
3269 <<constant tension parameter create/destroy function declarations>>=
3270 extern void *string_create_const_tension_param_t(char **param_strings);
3271 extern void destroy_const_tension_param_t(void *p);
3272 @
3273
3274 <<constant tension parameter create/destroy functions>>=
3275 const_tension_param_t *create_const_tension_param_t(double F)
3276 {
3277   const_tension_param_t *ret
3278     = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3279   assert(ret != NULL);
3280   ret->F = F;
3281   return ret;
3282 }
3283
3284 void *string_create_const_tension_param_t(char **param_strings)
3285 {
3286   return create_const_tension_param_t(atof(param_strings[0]));
3287 }
3288
3289 void destroy_const_tension_param_t(void *p)
3290 {
3291   if (p)
3292     free(p);
3293 }
3294
3295 @
3296
3297 <<constant tension model global declarations>>=
3298 extern int num_const_tension_params;
3299 extern char *const_tension_param_descriptions[];
3300 extern char const_tension_param_string[];
3301 @
3302
3303 <<constant tension model globals>>=
3304 int num_const_tension_params = 1;
3305 char *const_tension_param_descriptions[] = {"tension F, N"};
3306 char const_tension_param_string[] = "0";
3307 @
3308
3309 <<constant tension model getopt>>=
3310 {&const_tension_handler, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
3311 @
3312
3313 \subsection{Hooke}
3314 \label{sec.hooke}
3315
3316 <<hooke functions>>=
3317 <<hooke function>>
3318 <<hooke parameter create/destroy functions>>
3319 @
3320
3321 <<hooke tension model declarations>>=
3322 <<hooke tension function declaration>>
3323 <<hooke tension parameter create/destroy function declarations>>
3324 <<hooke tension model global declarations>>
3325
3326 @
3327
3328 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3329 The behavior of a series of springs $k_i$ in series is given by
3330 $$
3331   F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3332 $$
3333 For a simple proof, see Appendix \ref{app.math_hooke}.
3334
3335 <<hooke tension function declaration>>=
3336 extern double hooke_handler(double x, void *pdata);
3337 @
3338
3339 <<hooke function>>=
3340 double hooke_handler(double x, void *pdata)
3341 {
3342   tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3343   list_t *list = data->group;
3344   double k=0.0;
3345
3346   assert (x >= 0.0);
3347   list = head(list);
3348   assert(list != NULL); /* empty active group?! */
3349   while (list != NULL) {
3350     assert( ((hooke_param_t *)list->d)->k > 0 );
3351     k += 1.0/ ((hooke_param_t *)list->d)->k;
3352 #ifdef DEBUG
3353     fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3354             __FUNCTION__, ((hooke_param_t *)list->d)->k);
3355 #endif
3356     list = list->next;
3357   }
3358   k = 1.0 / k;
3359 #ifdef DEBUG
3360   fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
3361           __FUNCTION__, k, x, k*x);
3362 #endif
3363   return k*x;
3364 }
3365 @
3366
3367 <<hooke structure>>=
3368 typedef struct hooke_param_struct {
3369   double k; /* spring constant in N/m */
3370 } hooke_param_t;
3371 @
3372
3373 <<hooke tension parameter create/destroy function declarations>>=
3374 extern void *string_create_hooke_param_t(char **param_strings);
3375 extern void destroy_hooke_param_t(void *p);
3376 @
3377
3378 <<hooke parameter create/destroy functions>>=
3379 hooke_param_t *create_hooke_param_t(double k)
3380 {
3381   hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3382   assert(ret != NULL);
3383   ret->k = k;
3384   return ret;
3385 }
3386
3387 void *string_create_hooke_param_t(char **param_strings)
3388 {
3389   return create_hooke_param_t(atof(param_strings[0]));
3390 }
3391
3392 void destroy_hooke_param_t(void *p)
3393 {
3394   if (p)
3395     free(p);
3396 }
3397 @
3398
3399 <<hooke tension model global declarations>>=
3400 extern int num_hooke_params;
3401 extern char *hooke_param_descriptions[];
3402 extern char hooke_param_string[];
3403 @
3404
3405 <<hooke tension model globals>>=
3406 int num_hooke_params = 1;
3407 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3408 char hooke_param_string[]="0.05";
3409 @
3410
3411 <<hooke tension model getopt>>=
3412 {hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t}
3413 @
3414
3415 \subsection{Worm-like chain}
3416 \label{sec.wlc}
3417
3418 We can model several unfolded domains with a single worm-like chain.
3419 <<worm-like chain functions>>=
3420 <<worm-like chain function>>
3421 <<worm-like chain handler>>
3422 <<worm-like chain create/destroy functions>>
3423 @
3424
3425 <<worm-like chain tension model declarations>>=
3426 <<worm-like chain handler declaration>>
3427 <<worm-like chain create/destroy function declarations>>
3428 <<worm-like chain tension model global declarations>>
3429
3430 @
3431
3432 The combination of all unfolded domains is modeled as a worm like chain,
3433 with the force $F_{\text{WLC}}$ approximately given by
3434 $$
3435   F_{\text{WLC}} \approx \frac{k_B T}{p}
3436                \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3437 $$
3438 where $x$ is the distance between the fixed ends,
3439 $k_B$ is Boltzmann's constant,
3440 $T$ is the absolute temperature,
3441 $p$ is the persistence length, and
3442 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3443 <<worm-like chain function>>=
3444 double wlc(double x, double T, double p, double L)
3445 {/* N             m         K         m         m */
3446   double xL=0.0;               /* uses global k_B */
3447   assert(x >= 0);
3448   assert(T > 0); assert(p > 0); assert(L > 0);
3449   if (x >= L) return HUGE_VAL;
3450   xL = x/L; /* unitless */
3451   //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3452   //       k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3453   return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3454 }
3455 @
3456 This model assumes that each unfolded domain has the same persistence length.
3457
3458 <<worm-like chain handler declaration>>=
3459 extern double wlc_handler(double x, void *pdata);
3460 @
3461
3462 <<worm-like chain handler>>=
3463 double wlc_handler(double x, void *pdata)
3464 {
3465   tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3466   list_t *list = data->group;
3467   double p, L=0.0;
3468
3469   list = head(list);
3470   assert(list != NULL); /* empty active group?! */
3471   p = ((wlc_param_t *)list->d)->p;
3472   while (list != NULL) {
3473     /* enforce consistent p */
3474     assert( ((wlc_param_t *)list->d)->p == p);
3475     L += ((wlc_param_t *)list->d)->L;
3476 #ifdef DEBUG
3477     fprintf(stderr, "%s: WLC section %g m long in series\n",
3478             __FUNCTION__, ((wlc_param_t *)list->d)->L);
3479 #endif
3480     list = list->next;
3481   }
3482 #ifdef DEBUG
3483   fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3484           __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3485 #endif
3486   return wlc(x, data->env->T, p, L);
3487 }
3488 @
3489
3490 <<worm-like chain structure>>=
3491 typedef struct wlc_param_struct {
3492   double p;      /* WLC persistence length (m) */
3493   double L;      /* WLC contour length (m)     */
3494 } wlc_param_t;
3495 @
3496
3497 <<worm-like chain create/destroy function declarations>>=
3498 extern void *string_create_wlc_param_t(char **param_strings);
3499 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3500 @
3501
3502 <<worm-like chain create/destroy functions>>=
3503 wlc_param_t *create_wlc_param_t(double p, double L)
3504 {
3505   wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3506   assert(ret != NULL);
3507   ret->p = p;
3508   ret->L = L;
3509   return ret;
3510 }
3511
3512 void *string_create_wlc_param_t(char **param_strings)
3513 {
3514   return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3515 }
3516
3517 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3518 {
3519   if (p)
3520     free(p);
3521 }
3522 @
3523
3524 We define [[wlc_param_string]] as a global array because hard-coding the values into the tension model getopt structure gives a read-only version.
3525 TODO (now we copy the string before parsing, but still do this for clarity).
3526 For example,
3527 <<valid string write code>>=
3528 char string[] = "abc";
3529 string[1] = 'x';
3530 @ works, but
3531 <<valid string write code>>=
3532 char *string = "abc";
3533 string[1] = 'x';
3534 @ segfaults.
3535
3536 <<worm-like chain tension model global declarations>>=
3537 extern int num_wlc_params;
3538 extern char *wlc_param_descriptions[];
3539 extern char wlc_param_string[];
3540 @
3541 <<worm-like chain tension model globals>>=
3542 int num_wlc_params = 2;
3543 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3544 char wlc_param_string[]="0.39e-9,28.4e-9";
3545 @
3546
3547 <<worm-like chain tension model getopt>>=
3548 {wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
3549 @
3550 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3551
3552 \subsection{Freely-jointed chain}
3553 \label{sec.fjc}
3554
3555 <<freely-jointed chain functions>>=
3556 <<freely-jointed chain function>>
3557 <<freely-jointed chain handler>>
3558 <<freely-jointed chain create/destroy functions>>
3559 @
3560
3561 <<freely-jointed chain tension model declarations>>=
3562 <<freely-jointed chain handler declaration>>
3563 <<freely-jointed chain create/destroy function declarations>>
3564 <<freely-jointed chain tension model global declarations>>
3565
3566 @
3567
3568 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3569 The tension of a chain of $N$ such links is given by
3570 $$
3571   F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3572 $$
3573 where $L=Nl$ is the total length of the chain, and $\mathcal{L}(\alpha) \equiv \coth{\alpha} - \frac{1}{\alpha}$ is the Langevin function\citep{hatfield99}.
3574 <<freely-jointed chain function>>=
3575 double langevin(double x, void *params)
3576 {
3577   if (x==0) return 0;
3578   return 1.0/tanh(x) - 1/x;
3579 }
3580 <<one dimensional bracket declaration>>
3581 <<one dimensional solve declaration>>
3582 double inv_langevin(double x)
3583 {
3584   double lb=0.0, ub=1.0, ret;
3585   oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3586   ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3587   return ret;
3588 }
3589
3590 double fjc(double x, double T, double l, int N)
3591 {
3592   double L = l*(double)N;
3593   if (x == 0) return 0;
3594   //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3595   return k_B*T/l * inv_langevin(x/L);
3596 }
3597 @
3598
3599 <<freely-jointed chain handler declaration>>=
3600 extern double fjc_handler(double x, void *pdata);
3601 @
3602 <<freely-jointed chain handler>>=
3603 double fjc_handler(double x, void *pdata)
3604 {
3605   tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3606   list_t *list = data->group;
3607   double l;
3608   int N=0;
3609
3610   list = head(list);
3611   assert(list != NULL); /* empty active group?! */
3612   l = ((fjc_param_t *)list->d)->l;
3613   while (list != NULL) {
3614     /* enforce consistent l */
3615     assert( ((fjc_param_t *)list->d)->l == l);
3616     N += ((fjc_param_t *)list->d)->N;
3617 #ifdef DEBUG
3618     fprintf(stderr, "%s: FJC section %d links long in series\n",
3619             __FUNCTION__, ((fjc_param_t *)list->d)->N);
3620 #endif
3621     list = list->next;
3622   }
3623 #ifdef DEBUG
3624   fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3625           __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3626 #endif
3627   return fjc(x, data->env->T, l, N);
3628 }
3629 @
3630
3631 <<freely-jointed chain structure>>=
3632 typedef struct fjc_param_struct {
3633   double l;      /* FJC link length (m) */
3634   int N;         /* FJC number of links */
3635 } fjc_param_t;
3636 @
3637
3638 <<freely-jointed chain create/destroy function declarations>>=
3639 extern void *string_create_fjc_param_t(char **param_strings);
3640 extern void destroy_fjc_param_t(void *p);
3641 @
3642
3643 <<freely-jointed chain create/destroy functions>>=
3644 fjc_param_t *create_fjc_param_t(double l, double N)
3645 {
3646   fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3647   assert(ret != NULL);
3648   ret->l = l;
3649   ret->N = N;
3650   return ret;
3651 }
3652
3653 void *string_create_fjc_param_t(char **param_strings)
3654 {
3655   return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3656 }
3657
3658 void destroy_fjc_param_t(void *p)
3659 {
3660   if (p)
3661     free(p);
3662 }
3663 @
3664
3665 <<freely-jointed chain tension model global declarations>>=
3666 extern int num_fjc_params;
3667 extern char *fjc_param_descriptions[];
3668 extern char fjc_param_string[];
3669 @
3670
3671 <<freely-jointed chain tension model globals>>=
3672 int num_fjc_params = 2;
3673 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3674 char fjc_param_string[]="0.5e-9,200";
3675 @
3676
3677 <<freely-jointed chain tension model getopt>>=
3678 {fjc_handler, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
3679 @
3680
3681 \subsection{Tension model implementation}
3682
3683 <<tension-model.c>>=
3684 //#define DEBUG
3685 <<license comment>>
3686 <<tension model includes>>
3687 #include "tension_model.h"
3688 <<tension model internal definitions>>
3689 <<tension model globals>>
3690 <<tension model functions>>
3691 @
3692
3693 <<tension model includes>>=
3694 #include <assert.h> /* assert()                */
3695 #include <stdlib.h> /* NULL                    */
3696 #include <math.h>   /* HUGE_VAL, sqrt(), exp() */
3697 #include <stdio.h>  /* fprintf(), stdout       */
3698 #include "global.h"
3699 #include "list.h"
3700 #include "tension_balance.h" /* oneD_*() */
3701 @
3702
3703 <<tension model internal definitions>>=
3704 <<constant tension structure>>
3705 <<hooke structure>>
3706 <<worm-like chain structure>>
3707 <<freely-jointed chain structure>>
3708 @
3709
3710 <<tension model globals>>=
3711 <<tension handler array global>>
3712 <<constant tension model globals>>
3713 <<hooke tension model globals>>
3714 <<worm-like chain tension model globals>>
3715 <<freely-jointed chain tension model globals>>
3716 @
3717
3718 <<tension model functions>>=
3719 <<constant tension functions>>
3720 <<hooke functions>>
3721 <<worm-like chain functions>>
3722 <<freely-jointed chain functions>>
3723 @
3724
3725
3726 \subsection{Utilities}
3727
3728 The tension models can get complicated, and users may want to reassure
3729 themselves that this portion of the input physics are functioning properly. The
3730 stand-alone program [[tension_model_utils]] demonstrates the tension model
3731 interface without the overhead of having to understand a full simulation.
3732 [[tension_model_utils]] takes command line model arguments like the full
3733 simulation, and prints $F(x)$ data to the screen for the selected model over a
3734 range of $x$.
3735
3736 <<tension-model-utils.c>>=
3737 <<license comment>>
3738 <<tension model utility includes>>
3739 <<tension model utility definitions>>
3740 <<tension model utility getopt functions>>
3741 <<setup>>
3742 <<tension model utility main>>
3743 @
3744
3745 <<tension model utility main>>=
3746 int main(int argc, char **argv)
3747 {
3748   <<tension model getopt array definition>>
3749   tension_model_getopt_t *model = NULL;
3750   void *params;
3751   environment_t env;
3752   unsigned int flags;
3753   one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3754   tension_handler_data_t tdata;
3755   int num_param_args; /* for INIT_MODEL() */
3756   char **param_args;  /* for INIT_MODEL() */
3757   double Fmax,Xmax;
3758   get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
3759               &Fmax, &Xmax, &flags);
3760   setup();
3761   if (flags & VFLAG) {
3762     printf("#initializing model %s with parameters %s\n", model->name, model->params);
3763   }
3764   INIT_MODEL("utility", model, model->params, params);
3765   tdata.group = NULL;
3766   push(&tdata.group, params, 1);
3767   tdata.env = &env;
3768   tdata.persist = NULL;
3769   if (model->handler == NULL) {
3770     printf("No tension function for model %s\n", model->name);
3771     exit(0);
3772   }
3773   {
3774     int i,N=200;
3775     double x=0, F=0;
3776     printf("#Distance (m)\tForce (N)\n");
3777     for (i=0; i<=N; i++) {
3778       x = Xmax*i/(double)N;
3779       F = (*model->handler)(x, &tdata);
3780       if (F < 0 || F > Fmax) break;
3781       printf("%g\t%g\n", x, F);
3782     }
3783     if (flags & VFLAG && i <= N) {
3784       /* explain exit condition */
3785       if (F < 0)
3786         printf("Impossible force %g\n", F);
3787       else if (F > Fmax)
3788         printf("Reached large force limit %g > %g\n", F, Fmax);
3789     }
3790   }
3791   params = pop(&tdata.group);
3792   if (model->destructor)
3793     (*model->destructor)(params);
3794   return 0;
3795 }
3796 @
3797
3798 <<tension model utility includes>>=
3799 #include <stdlib.h>
3800 #include <stdio.h>
3801 #include <string.h> /* strlen() */
3802 #include <assert.h> /* assert() */
3803 #include "global.h"
3804 #include "parse.h"
3805 #include "list.h"
3806 #include "tension_model.h"
3807 @
3808
3809 <<tension model utility definitions>>=
3810 <<version definition>>
3811 #define VFLAG 1 /* verbose */
3812 <<string comparison definition and globals>>
3813 <<tension model getopt definitions>>
3814 <<initialize model definition>>
3815
3816 @
3817
3818 <<tension model utility getopt functions>>=
3819 <<index tension model>>
3820 <<tension model utility help>>
3821 <<tension model utility get options>>
3822 @
3823
3824 <<tension model utility help>>=
3825 void help(char *prog_name,
3826           environment_t *env,
3827           int n_tension_models, tension_model_getopt_t *tension_models,
3828           int tension_model, double Fmax, double Xmax)
3829 {
3830   int i, j;
3831   printf("usage: %s [options]\n", prog_name);
3832   printf("Version %s\n\n", VERSION);
3833   printf("A utility for understanding the available tension models\n\n");
3834   printf("Environmental options:\n");
3835   printf("-T T\tTemperature T (currently %g K)\n", env->T);
3836   printf("-C T\tYou can also set the temperature T in Celsius\n");
3837   printf("Model options:\n");
3838   printf("-m model\tFolded tension model (currently %s)\n",
3839          tension_models[tension_model].name);
3840   printf("-a args\tFolded tension model argument string (currently %s)\n",
3841          tension_models[tension_model].params);
3842   printf("F(x) is calculated for a range of x and printed\n");
3843   printf("For example:\n");
3844   printf("  #Distance (m)\tForce (N)\n");
3845   printf("  123.456\t7.89\n");
3846   printf("  ...\n");
3847   printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
3848   printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
3849   printf("-V\tChange output to verbose mode\n");
3850   printf("-h\tPrint this help and exit\n");
3851   printf("\n");
3852   printf("Tension models:\n");
3853   for (i=0; i<n_tension_models; i++) {
3854     printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3855     for (j=0; j < tension_models[i].num_params ; j++)
3856       printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3857     printf("\t\tdefaults: %s\n", tension_models[i].params);
3858   }
3859 }
3860 @
3861
3862 <<tension model utility get options>>=
3863 void get_options(int argc, char **argv, environment_t *env,
3864                  int n_tension_models, tension_model_getopt_t *tension_models,
3865                  tension_model_getopt_t **model, double *Fmax, double *Xmax,
3866                  unsigned int *flags)
3867 {
3868   char *prog_name = NULL;
3869   char c, options[] = "T:C:m:a:F:X:Vh";
3870   int tension_model=0, num_strings;
3871   extern char *optarg;
3872   extern int optind, optopt, opterr;
3873
3874   assert (argc > 0);
3875
3876   /* setup defaults */
3877
3878   prog_name = argv[0];
3879   env->T = 300.0;   /* K           */
3880   *Fmax = 1e5;
3881   *Xmax = 1e-6;
3882   *flags = 0;
3883   *model = tension_models;
3884
3885   while ((c=getopt(argc, argv, options)) != -1) {
3886     switch(c) {
3887     case 'T':  env->T   = atof(optarg);           break;
3888     case 'C':  env->T   = atof(optarg)+273.15;    break;
3889     case 'm':
3890       tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3891       *model = tension_models+tension_model;
3892       break;
3893     case 'a':
3894       tension_models[tension_model].params = optarg;
3895       break;
3896     case 'F': *Fmax = atof(optarg); break;
3897     case 'X': *Xmax = atof(optarg); break;
3898     case 'V': *flags |= VFLAG;      break;
3899     case '?':
3900       fprintf(stderr, "unrecognized option '%c'\n", optopt);
3901       /* fall through to default case */
3902     default:
3903       help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
3904       exit(1);
3905     }
3906   }
3907   return;
3908 }
3909 @
3910
3911
3912 \section{Unfolding rate models}
3913 \label{sec.k_models}
3914
3915 We have two main choices for determining $k$: Bell's model, which assumes the
3916 folded domains exist in a single `folded' state and make instantaneous,
3917 stochastic transitions over a free energy barrier; and Kramers' model, which
3918 treats the unfolding as a thermalized diffusion process.
3919 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3920 parameters into structures, with model specific functions for determining $k$.
3921
3922 <<k func definition>>=
3923 typedef double k_func_t(double F, environment_t *env, void *params);
3924
3925 @
3926 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3927 The interface is defined in [[k_model.h]], the implementation in [[k_model.c]], the unit testing in [[check_k_model.c]], and some other tidbits in [[k_model_utils.c]].
3928
3929 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3930 because \LaTeX\ doesn't like underscores outside math mode.
3931 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3932 You could use spaces instead (`\verb+ +'),
3933 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3934 which I don't like as much.
3935
3936 <<k-model.h>>=
3937 #ifndef K_MODEL_H
3938 #define K_MODEL_H
3939 <<license comment>>
3940 <<k func definition>>
3941 <<null k declarations>>
3942 <<const k declarations>>
3943 <<bell k declarations>>
3944 <<kbell k declarations>>
3945 <<kramers k declarations>>
3946 <<kramers integ k declarations>>
3947 #endif /* K_MODEL_H */
3948 @
3949
3950 <<k model module makefile lines>>=
3951 NW_SAWSIM_MODS += k_model
3952 CHECK_BINS += check_k_model
3953 check_k_model_MODS = parse string_eval
3954 @
3955 [[check_k_model]] also depends on the parse module.
3956
3957 <<k model utils makefile lines>>=
3958 K_MODEL_MODS = k_model parse string_eval
3959 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
3960         $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3961 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
3962         notangle -Rk-model-utils.c $< > $@
3963 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
3964         gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3965 clean_k_model_utils :
3966         rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
3967 @
3968
3969 \subsection{Null}
3970 \label{sec.null_k}
3971
3972 For domains that are never folded.
3973
3974 <<null k declarations>>=
3975 @
3976 <<null k functions>>=
3977 @
3978 <<null k globals>>=
3979 @
3980
3981 <<null k model getopt>>=
3982 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3983 @
3984
3985 \subsection{Constant rate model}
3986 \label{sec.k_const}
3987
3988 This is a pretty boring model, but it is usefull for testing the engine.
3989 $$
3990   k = k_0\;,
3991 $$
3992 where $k_0$ is the constant unfolding rate.
3993
3994 <<const k functions>>=
3995 <<const k function>>
3996 <<const k structure create/destroy functions>>
3997 @
3998
3999 <<const k declarations>>=
4000 <<const k function declaration>>
4001 <<const k structure create/destroy function declarations>>
4002 <<const k global declarations>>
4003
4004 @
4005
4006 <<const k structure>>=
4007 typedef struct const_k_param_struct {
4008   double knot;   /* transition rate k_0 (frac population per s) */
4009 } const_k_param_t;
4010 @
4011
4012 <<const k function declaration>>=
4013 double const_k(double F, environment_t *env, void *const_k_params);
4014 @
4015 <<const k function>>=
4016 double const_k(double F, environment_t *env, void *const_k_params)
4017 { /* Returns the rate constant k in frac pop per s. */
4018   const_k_param_t *p = (const_k_param_t *) const_k_params;
4019   assert(p != NULL);
4020   assert(p->knot > 0);
4021   return p->knot;
4022 }
4023 @
4024
4025 <<const k structure create/destroy function declarations>>=
4026 void *string_create_const_k_param_t(char **param_strings);
4027 void destroy_const_k_param_t(void *p);
4028 @
4029
4030 <<const k structure create/destroy functions>>=
4031 const_k_param_t *create_const_k_param_t(double knot)
4032 {
4033   const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4034   assert(ret != NULL);
4035   ret->knot = knot;
4036   return ret;
4037 }
4038
4039 void *string_create_const_k_param_t(char **param_strings)
4040 {
4041   return create_const_k_param_t(atof(param_strings[0]));
4042 }
4043
4044 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4045 {
4046   if (p)
4047     free(p);
4048 }
4049 @
4050
4051 <<const k global declarations>>=
4052 extern int num_const_k_params;
4053 extern char *const_k_param_descriptions[];
4054 extern char const_k_param_string[];
4055 @
4056
4057 <<const k globals>>=
4058 int num_const_k_params = 1;
4059 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4060 char const_k_param_string[]="1";
4061 @
4062
4063 <<const k model getopt>>=
4064 {"const", "constant rate transitions", &const_k, num_const_k_params, const_k_param_descriptions, (char *)const_k_param_string, &string_create_const_k_param_t, &destroy_const_k_param_t}
4065 @
4066
4067 \subsection{Bell's model}
4068 \label{sec.bell}
4069
4070 Quantitatively, Bell's model gives $k$ as
4071 $$
4072   k = k_0 \cdot e^{F dx / k_B T} \;,
4073 $$
4074 where $k_0$ is the unforced unfolding rate,
4075 $F$ is the applied unfolding force,
4076 $dx$ is the distance to the transition state, and
4077 $k_B T$ is the thermal energy\citep{schlierf06}.
4078
4079 <<bell k functions>>=
4080 <<bell k function>>
4081 <<bell k structure create/destroy functions>>
4082 @
4083
4084 <<bell k declarations>>=
4085 <<bell k function declaration>>
4086 <<bell k structure create/destroy function declarations>>
4087 <<bell k global declarations>>
4088
4089 @
4090
4091 <<bell k structure>>=
4092 typedef struct bell_param_struct {
4093   double knot;   /* transition rate k_0 (frac population per s) */
4094   double dx;     /* `distance to transition state' (nm) */
4095 } bell_param_t;
4096 @
4097
4098 <<bell k function declaration>>=
4099 double bell_k(double F, environment_t *env, void *bell_params);
4100 @
4101 <<bell k function>>=
4102 double bell_k(double F, environment_t *env, void *bell_params)
4103 { /* Returns the rate constant k in frac pop per s.
4104    * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4105    * uses global k_B in J/K */
4106   bell_param_t *p = (bell_param_t *) bell_params;
4107   assert(F >= 0); assert(env->T > 0);
4108   assert(p != NULL);
4109   assert(p->knot > 0); assert(p->dx >= 0);
4110 /*
4111   printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4112          F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4113          p->knot * exp(F*p->dx / (k_B*env->T) ));
4114   printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4115 */
4116   return p->knot * exp(F*p->dx / (k_B*env->T) );
4117 }
4118 @
4119
4120 <<bell k structure create/destroy function declarations>>=
4121 void *string_create_bell_param_t(char **param_strings);
4122 void destroy_bell_param_t(void *p);
4123 @
4124
4125 <<bell k structure create/destroy functions>>=
4126 bell_param_t *create_bell_param_t(double knot, double dx)
4127 {
4128   bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4129   assert(ret != NULL);
4130   ret->knot = knot;
4131   ret->dx = dx;
4132   return ret;
4133 }
4134
4135 void *string_create_bell_param_t(char **param_strings)
4136 {
4137   return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4138 }
4139
4140 void destroy_bell_param_t(void *p /* bell_param_t * */)
4141 {
4142   if (p)
4143     free(p);
4144 }
4145 @
4146
4147 <<bell k global declarations>>=
4148 extern int num_bell_params;
4149 extern char *bell_param_descriptions[];
4150 extern char bell_param_string[];
4151 @
4152
4153 <<bell k globals>>=
4154 int num_bell_params = 2;
4155 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4156 char bell_param_string[]="3.3e-4,0.25e-9";
4157 @
4158
4159 <<bell k model getopt>>=
4160 {"bell", "thermalized, two-state transitions", &bell_k, num_bell_params, bell_param_descriptions, (char *)bell_param_string, &string_create_bell_param_t, &destroy_bell_param_t}
4161 @
4162 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4163 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4164
4165
4166 \subsection{Linker stiffness corrected Bell model}
4167 \label{sec.kbell}
4168
4169 Walton et. al showed that the Bell model constant force approximation
4170 doesn't hold for biotin-streptavadin binding, and I extended their
4171 results to I27 for the 2009 Biophysical Society Annual
4172 Meeting\cite{walton08,king09}.  More details to follow when I get done
4173 with the conference...
4174
4175 We adjust Bell's model to
4176 $$
4177   k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4178 $$
4179 where $k_0$ is the unforced unfolding rate,
4180 $F$ is the applied unfolding force,
4181 $\kappa$ is the effective spring constant,
4182 $dx$ is the distance to the transition state, and
4183 $k_B T$ is the thermal energy\citep{schlierf06}.
4184
4185 <<kbell k functions>>=
4186 <<kbell k function>>
4187 <<kbell k structure create/destroy functions>>
4188 @
4189
4190 <<kbell k declarations>>=
4191 <<kbell k function declaration>>
4192 <<kbell k structure create/destroy function declarations>>
4193 <<kbell k global declarations>>
4194
4195 @
4196
4197 <<kbell k structure>>=
4198 typedef struct kbell_param_struct {
4199   double knot;   /* transition rate k_0 (frac population per s) */
4200   double dx;     /* `distance to transition state' (nm) */
4201 } kbell_param_t;
4202 @
4203
4204 <<kbell k function declaration>>=
4205 double kbell_k(double F, environment_t *env, void *kbell_params);
4206 @
4207 <<kbell k function>>=
4208 double kbell_k(double F, environment_t *env, void *kbell_params)
4209 { /* Returns the rate constant k in frac pop per s.
4210    * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4211    * uses global k_B in J/K */
4212   kbell_param_t *p = (kbell_param_t *) kbell_params;
4213   assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4214   assert(p != NULL);
4215   assert(p->knot > 0); assert(p->dx >= 0);
4216   return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4217 }
4218 @
4219
4220 <<kbell k structure create/destroy function declarations>>=
4221 void *string_create_kbell_param_t(char **param_strings);
4222 void destroy_kbell_param_t(void *p);
4223 @
4224
4225 <<kbell k structure create/destroy functions>>=
4226 kbell_param_t *create_kbell_param_t(double knot, double dx)
4227 {
4228   kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4229   assert(ret != NULL);
4230   ret->knot = knot;
4231   ret->dx = dx;
4232   return ret;
4233 }
4234
4235 void *string_create_kbell_param_t(char **param_strings)
4236 {
4237   return create_kbell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4238 }
4239
4240 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4241 {
4242   if (p)
4243     free(p);
4244 }
4245 @
4246
4247 <<kbell k global declarations>>=
4248 extern int num_kbell_params;
4249 extern char *kbell_param_descriptions[];
4250 extern char kbell_param_string[];
4251 @
4252
4253 <<kbell k globals>>=
4254 int num_kbell_params = 2;
4255 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4256 char kbell_param_string[]="3.3e-4,0.25e-9";
4257 @
4258
4259 <<kbell k model getopt>>=
4260 {"kbell", "thermalized, two-state transitions, corrected for effective linker stiffness", &kbell_k, num_kbell_params, kbell_param_descriptions, (char *)kbell_param_string, &string_create_kbell_param_t, &destroy_kbell_param_t}
4261 @
4262 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4263 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4264
4265
4266 \subsection{Kramer's model}
4267 \label{sec.kramers}
4268
4269 Kramer's model gives $k$ as
4270 %$$
4271 %  \frac{1}{k} = \frac{1}{D}
4272 %                \int_{x_\text{min}}^{x_\text{max}}
4273 %                     e^{\frac{-U_F(x)}{k_B T}}
4274 %                     \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4275 %                     dx,
4276 %$$
4277 %where $D$ is the diffusion constant,
4278 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4279 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4280 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4281 $$
4282   k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4283 $$
4284 where $D$ is the diffusion constant,
4285 $$
4286   l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4287 $$
4288 are length scales where
4289 $x_c(F)$ is the location of the bound state and
4290 $x_{ts}(F)$ is the location of the transition state,
4291 $E(F,x)$ is the free energy, and
4292 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4293 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4294  \citet{evans97} Eqn.~3).
4295
4296 <<kramers k functions>>=
4297 <<kramers k function>>
4298 <<kramers k structure create/destroy functions>>
4299 @
4300
4301 <<kramers k declarations>>=
4302 <<kramers k function declaration>>
4303 <<kramers k structure create/destroy function declarations>>
4304 <<kramers k global declarations>>
4305
4306 @
4307
4308 <<kramers k structure>>=
4309 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4310 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4311
4312 typedef struct kramers_param_struct {
4313   double D;                 /* diffusion rate (um/s)                 */
4314   char *xc;
4315   char *xts;
4316   char *ddEdxx;
4317   char *E;
4318   FILE *in;
4319   FILE *out;
4320   //kramers_x_func_t *xc;     /* function returning metastable x (nm)  */
4321   //kramers_x_func_t *xts;    /* function returning transition x (nm)  */
4322   //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4323   //kramers_E_func_t *E;      /* function returning E (J)              */
4324   //double *E_params;         /* parametrize the energy landscape     */
4325   //int n_E_params;           /* length of E_params array              */
4326 } kramers_param_t;
4327 @
4328
4329 <<kramers k function declaration>>=
4330 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4331 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4332 @
4333 <<kramers k function>>=
4334 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4335 {
4336   static char input[160]={0};
4337   static char output[80]={0};
4338   /* setup the environment */
4339   fprintf(in, "F = %g; T = %g\n", F, T);
4340   sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4341   string_eval(in, out, input, 80, output);
4342   return atof(output);
4343 }
4344
4345 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4346 {
4347   static char input[160]={0};
4348   static char output[80]={0};
4349   /* setup the environment */
4350   fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4351   sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4352   string_eval(in, out, input, 80, output);
4353   return atof(output);
4354 }
4355
4356 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4357 {
4358   kramers_param_t *p = (kramers_param_t *) kramers_params;
4359   return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4360 }
4361
4362 double kramers_k(double F, environment_t *env, void *kramers_params)
4363 { /* Returns the rate constant k in frac pop per s.
4364    * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4365    * uses global k_B in J/K */
4366   kramers_param_t *p = (kramers_param_t *) kramers_params;
4367   double kbT, xc, xts, lc, lts, Eb;
4368   assert(F >= 0); assert(env->T > 0);
4369   assert(p != NULL);
4370   assert(p->D > 0);
4371   assert(p->xc != NULL); assert(p->xts != NULL);
4372   assert(p->ddEdxx != NULL); assert(p->E != NULL);
4373   assert(p->in != NULL); assert(p->out != NULL);
4374   kbT = k_B*env->T;
4375   xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4376   if (xc == -1.0) return -1.0; /* error (Too much force) */
4377   xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4378   if (xc == -1.0) return -1.0; /* error (Too much force) */
4379   lc = sqrt(2.0*M_PI*kbT /
4380             kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4381   lts = sqrt(-2.0*M_PI*kbT/
4382             kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4383   Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4384      - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4385   //fprintf(stderr,"D = %g, lc = %g, lts = %g, Eb = %g, kbT = %g, k = %g\n", p->D, lc, lts, Eb, kbT, p->D/(lc*lts) * exp(-Eb/kbT));
4386   return p->D/(lc*lts) * exp(-Eb/kbT);
4387 }
4388 @
4389
4390 <<kramers k structure create/destroy function declarations>>=
4391 void *string_create_kramers_param_t(char **param_strings);
4392 void destroy_kramers_param_t(void *p);
4393 @
4394
4395 <<kramers k structure create/destroy functions>>=
4396 kramers_param_t *create_kramers_param_t(double D,
4397                                         char *xc_fn, char *xts_fn,
4398                                         char *E_fn, char *ddEdxx_fn,
4399                                         char *global_define_string)
4400 //                              kramers_x_func_t *xc, kramers_x_func_t *xts,
4401 //                            kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4402 //                            double *E_params, int n_E_params)
4403 {
4404   kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4405   assert(ret != NULL);
4406   ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4407   assert(ret->xc != NULL);
4408   ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4409   assert(ret->xts != NULL);
4410   ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4411   assert(ret->E != NULL);
4412   ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4413   assert(ret->ddEdxx != NULL);
4414   ret->D = D;
4415   strcpy(ret->xc, xc_fn);
4416   strcpy(ret->xts, xts_fn);
4417   strcpy(ret->E, E_fn);
4418   strcpy(ret->ddEdxx, ddEdxx_fn);
4419   //ret->E_params = E_params;
4420   //ret->n_E_params = n_E_params;
4421   string_eval_setup(&ret->in, &ret->out);
4422   fprintf(ret->in, "kB = %g\n", k_B);
4423   fprintf(ret->in, "%s\n", global_define_string);
4424   return ret;
4425 }
4426
4427 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4428 void *string_create_kramers_param_t(char **param_strings)
4429 {
4430   return create_kramers_param_t(atof(param_strings[0]),
4431                                 param_strings[2],
4432                                 param_strings[3],
4433                                 param_strings[4],
4434                                 param_strings[5],
4435                                 param_strings[1]);
4436 }
4437
4438 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4439 {
4440   kramers_param_t *pk = (kramers_param_t *)p;
4441   if (pk) {
4442     string_eval_teardown(&pk->in, &pk->out);
4443     //if (pk->E_params)
4444     //  free(pk->E_params);
4445     free(pk);
4446   }
4447 }
4448 @
4449
4450 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4451 Schlierf and Rief used a cubic-spline interpolation routine and the double integral form of Kramers' theory, so we get to pick an actual function to approximate the energy landscape.
4452 We follow \citet{shillcock98} and use
4453 $$
4454   E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4455                          \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4456                         -F x
4457 $$
4458 where TODO, variable meanings.%\citep{shillcock98}.
4459 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4460 \begin{align}
4461   E'(F,E_b,x,x_s) &=\frac{27 E_b}{4 x_s}\frac{x}{x_s}\p({4-9\frac{x}{x_s}})-F\\
4462   E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4463 \end{align}
4464 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4465 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4466 The roots are therefor at
4467 \begin{align}
4468   x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4469         &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4470         &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4471 \end{align}
4472
4473 <<kramers k global declarations>>=
4474 extern int num_kramers_params;
4475 extern char *kramers_param_descriptions[];
4476 extern char kramers_param_string[];
4477 @
4478
4479 <<kramers k globals>>=
4480 int num_kramers_params = 6;
4481 char *kramers_param_descriptions[] = {"Diffusion constant D, m^2/s", "constant parameter declarations", "bound position xc, m", "transition state xts, m","energy E, J", "d^2E/dx^2, J/m^2"};
4482 char kramers_param_string[]="0.2e-12,{Eb = 20*300*kB;xs = 1.5e-9},{(F*xs > 3*Eb) and (-1) or (2*xs/9.0*(1 - (1 - F*xs/(3*Eb))**0.5))},{2*xs/9.0*(1 + (1 - F*xs/(3*Eb))**0.5)},{27.0/4 * Eb * (x/xs)**2 * (2-3*x/xs) - F*x},{27.0*Eb/(2*xs**2) * (2-9*x/xs)}";
4483 @
4484 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4485 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4486 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4487 Returning an $x_c$ position of $-1\U{m}$ lets the simulation know that the rate is poorly defined for that force, and the timestepping algorithm in Section \ref{app.adaptive_dt} reduces it's timestep appropriately.
4488 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4489 It works so long as [[val_1]] is not [[false]].
4490
4491 <<kramers k model getopt>>=
4492 {"kramers", "thermalized, diffusion-limited transitions (function parameters are parsed by Python.  For distances, the environment variables force 'F' and temperature 'T' are defined, as well as the Boltzmann constant 'kB'.  For energies, the position 'x' is also defined.  Functions must evaluate to a float representing their output and produce no output on stdout.", &kramers_k, num_kramers_params, kramers_param_descriptions, (char *)kramers_param_string, &string_create_kramers_param_t, &destroy_kramers_param_t}
4493 @
4494
4495 \subsection{Kramer's model (natural cubic splines)}
4496 \label{sec.kramers_integ}
4497
4498 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4499 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4500 \citet{schlierf06} Eqn.~2)
4501 $$
4502   \frac{1}{k} = \frac{1}{D}
4503                 \int_{x_\text{min}}^{x_\text{max}}
4504                      e^{\frac{U_F(x)}{k_B T}}
4505                      \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4506                      dx,
4507 $$
4508 where $D$ is the diffusion constant,
4509 $U_F(x)$ is the free energy along the unfolding coordinate, and
4510 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4511  before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4512
4513 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4514 interpolating between them with cubic splines.
4515 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4516
4517 <<kramers integ k functions>>=
4518 <<spline functions>>
4519 <<kramers integ A k functions>>
4520 <<kramers integ B k functions>>
4521 <<kramers integ k function>>
4522 <<kramers integ k structure create/destroy functions>>
4523 @
4524
4525 <<kramers integ k declarations>>=
4526 <<kramers integ k function declaration>>
4527 <<kramers integ k structure create/destroy function declarations>>
4528 <<kramers integ k global declarations>>
4529
4530 @
4531
4532 <<kramers integ k structure>>=
4533 typedef double func_t(double x, void *params);
4534 typedef struct kramers_integ_param_struct {
4535   double D;            /* diffusion rate (um/s)                 */
4536   double x_min;        /* integration bounds                    */
4537   double x_max;
4538   func_t *E;           /* function returning E (J)              */
4539   void *E_params;      /* parametrize the energy landscape     */
4540   destroy_data_func_t *destroy_E_params;
4541
4542   double F;            /* for passing into GSL evaluations      */
4543   environment_t *env;
4544 } kramers_integ_param_t;
4545 @
4546
4547 <<spline param structure>>=
4548 typedef struct spline_param_struct {
4549   int n_knots;            /* natural cubic spline knots for U(x)   */
4550   double *x;
4551   double *y;
4552   gsl_spline *spline;    /* GSL spline parameters               */
4553   gsl_interp_accel *acc; /* GSL spline acceleration data        */
4554 } spline_param_t;
4555 @
4556
4557 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4558 $$
4559   \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4560 $$
4561 <<kramers integ A k functions>>=
4562 double kramers_integ_k_integrandA(double x, void *params)
4563 {
4564   kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4565   double U = (*p->E)(x, p->E_params);
4566   double Fx = p->F * x;
4567   double kbT = k_B * p->env->T;
4568   //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4569   return exp(-(U-Fx)/kbT);
4570 }
4571 double kramers_integ_k_integralA(double x, void *params)
4572 {
4573   kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4574   gsl_function f;
4575   double result, abserr;
4576   size_t neval; /* for qng */
4577   static gsl_integration_workspace *w = NULL;
4578   if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4579   f.function = &kramers_integ_k_integrandA;
4580   f.params = params;
4581   //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4582   assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4583   //fprintf(stderr, "integralA = %g\n", result);
4584   return result;
4585 }
4586 @
4587
4588 $$
4589   \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4590                 \int_{x_\text{min}}^{x_\text{max}}
4591                      e^{\frac{U_F(x)}{k_B T}}
4592                      \text{Integral}_A(x)
4593                      dx,
4594 $$
4595 <<kramers integ B k functions>>=
4596 double kramers_integ_k_integrandB(double x, void *params)
4597 {
4598   kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4599   double U = (*p->E)(x, p->E_params);
4600   double Fx = p->F * x;
4601   double kbT = k_B * p->env->T;
4602   double intA = kramers_integ_k_integralA(x, params);
4603   //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4604   return exp((U-Fx)/kbT)*intA;
4605 }
4606 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4607 {
4608   kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4609   gsl_function f;
4610   double result, abserr;
4611   size_t neval; /* for qng */
4612   static gsl_integration_workspace *w = NULL;
4613   if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4614   f.function = &kramers_integ_k_integrandB;
4615   f.params = params;
4616   p->F = F;
4617   p->env = env;
4618   //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4619   assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4620   //fprintf(stderr, "integralB = %g\n", result);
4621   return result;
4622 }
4623 @
4624
4625 With the integrals taken care of, evaluating $k$ is simply
4626 $$
4627   k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4628 $$
4629 <<kramers integ k function declaration>>=
4630 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4631 @
4632 <<kramers integ k function>>=
4633 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4634 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4635   kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4636   return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4637 }
4638 @
4639
4640 <<kramers integ k structure create/destroy function declarations>>=
4641 void *string_create_kramers_integ_param_t(char **param_strings);
4642 void destroy_kramers_integ_param_t(void *p);
4643 @
4644
4645 <<kramers integ k structure create/destroy functions>>=
4646 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4647                               double xmin, double xmax, func_t *E,
4648                               void *E_params,
4649                               destroy_data_func_t *destroy_E_params)
4650 {
4651   kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4652   assert(ret != NULL);
4653   ret->D = D;
4654   ret->x_min = xmin;
4655   ret->x_max = xmax;
4656   ret->E = E;
4657   ret->E_params = E_params;
4658   ret->destroy_E_params = destroy_E_params;
4659   return ret;
4660 }
4661
4662 void *string_create_kramers_integ_param_t(char **param_strings)
4663 {
4664   //printf("create kramers integ.  D: %s, knots: %s, x in [%s, %s]\n",
4665   //       param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4666   void *E_params = string_create_spline_param_t(param_strings+1);
4667   return create_kramers_integ_param_t(atof(param_strings[0]),
4668                                       atof(param_strings[2]),
4669                                       atof(param_strings[3]),
4670                                       &spline_eval, E_params,
4671                                       destroy_spline_param_t);
4672 }
4673
4674 void destroy_kramers_integ_param_t(void *params)
4675 {
4676   kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4677   if (p) {
4678     if (p->E_params)
4679       (*p->destroy_E_params)(p->E_params);
4680     free(p);
4681   }
4682 }
4683 @
4684
4685 Finally we have the GSL spline wrappers:
4686 <<spline functions>>=
4687 <<spline function>>
4688 <<create/destroy spline>>
4689 @
4690
4691 <<spline function>>=
4692 double spline_eval(double x, void *spline_params)
4693 {
4694   spline_param_t *p = (spline_param_t *)(spline_params);
4695   return gsl_spline_eval(p->spline, x, p->acc);
4696 }
4697 @
4698
4699 <<create/destroy spline>>=
4700 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4701 {
4702   spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4703   assert(ret != NULL);
4704   ret->n_knots = n_knots;
4705   ret->x = x;
4706   ret->y = y;
4707   ret->acc = gsl_interp_accel_alloc();
4708   ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4709   gsl_spline_init(ret->spline, x, y, n_knots);
4710   return ret;
4711 }
4712
4713 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4714 void *string_create_spline_param_t(char **param_strings)
4715 {
4716   char **knot_coords;
4717   int i, num_knots;
4718   double *x=NULL, *y=NULL;
4719   /* break into ordered pairs */
4720   parse_list_string(param_strings[0], SEP, '(', ')',
4721                     &num_knots, &knot_coords);
4722   x = (double *)malloc(sizeof(double)*num_knots);
4723   assert(x != NULL);
4724   y = (double *)malloc(sizeof(double)*num_knots);
4725   assert(x != NULL);
4726   for (i=0; i<num_knots; i++) {
4727     if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4728       fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4729       exit(1);
4730     }
4731     //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4732   }
4733   free_parsed_list(num_knots, knot_coords);
4734   return create_spline_param_t(num_knots, x, y);
4735 }
4736
4737 void destroy_spline_param_t(void *params)
4738 {
4739   spline_param_t *p = (spline_param_t *)params;
4740   if (p) {
4741     if (p->spline)
4742       gsl_spline_free(p->spline);
4743     if (p->acc)
4744       gsl_interp_accel_free(p->acc);
4745     if (p->x)
4746       free(p->x);
4747     if (p->y)
4748       free(p->y);
4749     free(p);
4750   }
4751 }
4752 @
4753
4754 <<kramers integ k global declarations>>=
4755 extern int num_kramers_integ_params;
4756 extern char *kramers_integ_param_descriptions[];
4757 extern char kramers_integ_param_string[];
4758 @
4759
4760 <<kramers integ k globals>>=
4761 int num_kramers_integ_params = 4;
4762 char *kramers_integ_param_descriptions[] = {
4763                            "diffusion D, m^2 / s",
4764                            "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4765                            "starting integration bound x_min, m",
4766                            "ending integration bound x_max, m"};
4767 //char kramers_integ_param_string[]="0.2e-12,{(3.5e-9,-5*k_B*300),(4e-9,-20*k_B*300),(4.25e-9,-8*k_B*300),(5e-9,0),(5.5e-9,-10*k_B*300)},3e-9,6e-9";
4768 //char kramers_integ_param_string[]="0.2e-12,{(3.5e-9,-2.1e-20),(4e-9,-8.3e-20),(4.25e-9,-3.3e-20),(5e-9,0),(5.5e-9,-4.1e-20)},3e-9,6e-9";
4769 char kramers_integ_param_string[]="0.2e-12,{(2e-9,0),(2.1e-9,-3.7e-20),(2.2e-9,-4.4e-20),(2.35e-9,-2.5e-20),(2.41e-9,-7e-21),(2.45e-9,8e-22),(2.51e-9,6.9e-21),(2.64e-9,1.39e-20),(2.9e-9,2.55e-20),(3.52e-9,2.9e-20),(3.7e-9,1.45e-20),(4e-9,-1.7e-20),(6e-9,-2e-20),(9e-9,-3e-20),(1e-8,-3e-20)},2e-9,6e-9";
4770 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4771  * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4772  * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4773 @
4774
4775 <<kramers integ k model getopt>>=
4776 {"kramers_integ", "thermalized, diffusion-limited transitions", &kramers_integ_k, num_kramers_integ_params, kramers_integ_param_descriptions, (char *)kramers_integ_param_string, &string_create_kramers_integ_param_t, &destroy_kramers_integ_param_t}
4777 @
4778 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4779 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4780
4781 \subsection{Unfolding model implementation}
4782
4783 <<k-model.c>>=
4784 <<license comment>>
4785 <<k model includes>>
4786 #include "k_model.h"
4787 <<k model internal definitions>>
4788 <<k model globals>>
4789 <<k model functions>>
4790 @
4791
4792 <<k model includes>>=
4793 #include <assert.h> /* assert()                */
4794 #include <stdlib.h> /* NULL, malloc()          */
4795 #include <math.h>   /* HUGE_VAL, sqrt(), exp() */
4796 #include <stdio.h>  /* fprintf(), stdout       */
4797 #include <string.h> /* strlen(), strcpy()      */
4798 #include <gsl/gsl_integration.h>
4799 #include <gsl/gsl_spline.h>
4800 #include "global.h"
4801 #include "parse.h"
4802 @
4803
4804 <<k model internal definitions>>=
4805 <<const k structure>>
4806 <<bell k structure>>
4807 <<kbell k structure>>
4808 <<kramers k structure>>
4809 <<kramers integ k structure>>
4810 <<spline param structure>>
4811 @
4812
4813 <<k model globals>>=
4814 <<null k globals>>
4815 <<const k globals>>
4816 <<bell k globals>>
4817 <<kbell k globals>>
4818 <<kramers k globals>>
4819 <<kramers integ k globals>>
4820 @
4821
4822 <<k model functions>>=
4823 <<null k functions>>
4824 <<const k functions>>
4825 <<bell k functions>>
4826 <<kbell k functions>>
4827 <<kramers k functions>>
4828 <<kramers integ k functions>>
4829 @
4830
4831 \subsection{Unfolding model unit tests}
4832
4833 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4834 <<check-k-model.c>>=
4835 <<k model check includes>>
4836 <<check relative error>>
4837 <<model globals>>
4838 <<k model test suite>>
4839 <<main check program>>
4840 @
4841
4842 <<k model check includes>>=
4843 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4844 #include <stdio.h>  /* sprintf()                             */
4845 #include <assert.h> /* assert()                              */
4846 #include <math.h>   /* exp()                                 */
4847 <<check includes>>
4848 #include "global.h"
4849 #include "k_model.h"
4850 @
4851
4852 <<k model test suite>>=
4853 <<const k tests>>
4854 <<bell k tests>>
4855 <<k model suite definition>>
4856 @
4857
4858 <<k model suite definition>>=
4859 Suite *test_suite (void)
4860 {
4861   Suite *s = suite_create ("k model");
4862   <<const k test case defs>>
4863   <<bell k test case defs>>
4864
4865   <<const k test case adds>>
4866   <<bell k test case adds>>
4867   return s;
4868 }
4869 @
4870
4871 \subsubsection{Constant}
4872
4873 <<const k test case defs>>=
4874 TCase *tc_const_k = tcase_create("Constant k");
4875 @
4876
4877 <<const k test case adds>>=
4878 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4879 tcase_add_test(tc_const_k, test_const_k_over_range);
4880 suite_add_tcase(s, tc_const_k);
4881 @
4882
4883
4884 <<const k tests>>=
4885 START_TEST(test_const_k_create_destroy)
4886 {
4887   char *knot[] = {"1","2","3","4","5","6"};
4888   char *params[] = {knot[0]};
4889   void *p = NULL;
4890   int i;
4891   for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4892     params[0] = knot[i];
4893     p = string_create_const_k_param_t(params);
4894     destroy_const_k_param_t(p);
4895   }
4896 }
4897 END_TEST
4898
4899 START_TEST(test_const_k_over_range)
4900 {
4901   double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4902   char *knot[] = {"1","2","3","4","5","6"};
4903   char *params[] = {knot[0]};
4904   void *p = NULL;
4905   environment_t env;
4906   char param_string[80];
4907   int i;
4908   env.T = T;
4909   for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4910     params[0] = knot[i];
4911     p = string_create_const_k_param_t(params);
4912     for ( F=Fm; F<FM; F+=dF ) {
4913       fail_unless(const_k(F, &env, p)==atof(knot[i]),
4914           "const_k(%g, %g, &{%s}) = %g != %s",
4915           F, T, knot[i], const_k(F, &env, p), knot[i]);
4916     }
4917     destroy_const_k_param_t(p);
4918   }
4919 }
4920 END_TEST
4921 @
4922
4923 \subsubsection{Bell}
4924
4925 <<bell k test case defs>>=
4926 TCase *tc_bell = tcase_create("Bell's k");
4927 @
4928
4929 <<bell k test case adds>>=
4930 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4931 tcase_add_test(tc_bell, test_bell_k_at_zero);
4932 tcase_add_test(tc_bell, test_bell_k_at_one);
4933 suite_add_tcase(s, tc_bell);
4934 @
4935
4936 <<bell k tests>>=
4937 START_TEST(test_bell_k_create_destroy)
4938 {
4939   char *knot[] = {"1","2","3","4","5","6"};
4940   char *dx="1";
4941   char *params[] = {knot[0], dx};
4942   void *p = NULL;
4943   int i;
4944   for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4945     params[0] = knot[i];
4946     p = string_create_bell_param_t(params);
4947     destroy_bell_param_t(p);
4948   }
4949 }
4950 END_TEST
4951
4952 START_TEST(test_bell_k_at_zero)
4953 {
4954   double F=0.0, T=300.0;
4955   char *knot[] = {"1","2","3","4","5","6"};
4956   char *dx="1";
4957   char *params[] = {knot[0], dx};
4958   void *p = NULL;
4959   environment_t env;
4960   char param_string[80];
4961   int i;
4962   env.T = T;
4963   for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4964     params[0] = knot[i];
4965     p = string_create_bell_param_t(params);
4966     fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4967                 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4968                 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4969     destroy_bell_param_t(p);
4970   }
4971 }
4972 END_TEST
4973
4974 START_TEST(test_bell_k_at_one)
4975 {
4976   double T=300.0;
4977   char *knot[] = {"1","2","3","4","5","6"};
4978   char *dx="1";
4979   char *params[] = {knot[0], dx};
4980   double F= k_B*T/atof(dx);
4981   void *p = NULL;
4982   environment_t env;
4983   char param_string[80];
4984   int i;
4985   env.T = T;
4986   for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4987     params[0] = knot[i];
4988     p = string_create_bell_param_t(params);
4989     CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4990     //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4991     //            "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4992     //            F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4993     destroy_bell_param_t(p);
4994   }
4995 }
4996 END_TEST
4997 @
4998
4999 <<kramers k tests>>=
5000 @
5001
5002 <<kramers k test case def>>=
5003 @
5004
5005 <<kramers k test case add>>=
5006 @
5007
5008 <<k model function tests>>=
5009 <<check relative error>>
5010
5011 START_TEST(test_something)
5012 {
5013   double k=5, x=3, last_x=2.0, F;
5014   one_dim_fn_t *handlers[] = {&hooke};
5015   void *data[] =               {&k};
5016   double xi[] =                {0.0};
5017   int active[] =               {1};
5018   int new_active_groups = 1;
5019   F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5020   fail_unless(F = k*x, NULL);
5021 }
5022 END_TEST
5023 @
5024
5025 \subsection{Utilities}
5026
5027 The unfolding models can get complicated, and are sometimes hard to explain to others.
5028 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5029 the overhead of having to understand a full simulation.
5030 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5031 or special mode, where the operation depends on the specific model selected.
5032
5033 <<k-model-utils.c>>=
5034 <<license comment>>
5035 <<k model utility includes>>
5036 <<k model utility definitions>>
5037 <<k model utility getopt functions>>
5038 <<k model utility multi model E>>
5039 <<k model utility main>>
5040 @
5041
5042 <<k model utility main>>=
5043 int main(int argc, char **argv)
5044 {
5045   <<k model getopt array definition>>
5046   k_model_getopt_t *model = NULL;
5047   void *params;
5048   enum mode_t mode;
5049   environment_t env;
5050   unsigned int flags;
5051   int num_param_args; /* for INIT_MODEL() */
5052   char **param_args;  /* for INIT_MODEL() */
5053   double Fmax;
5054   double special_xmin;
5055   double special_xmax;
5056   get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5057               &Fmax, &special_xmin, &special_xmax, &flags);
5058   if (flags & VFLAG) {
5059     printf("#initializing model %s with parameters %s\n", model->name, model->params);
5060   }
5061   INIT_MODEL("utility", model, model->params, params);
5062   switch (mode) {
5063     case M_K_OF_F :
5064       if (model->k == NULL) {
5065         printf("No k function for model %s\n", model->name);
5066         exit(0);
5067       }
5068       if (Fmax <= 0) {
5069         printf("Fmax = %g <= 0\n", Fmax);
5070         exit(1);
5071       }
5072       {
5073         double F, k=1.0;
5074         int i, N=200;
5075         printf("#F (N)\tk (%% pop. per s)\n");
5076         for (i=0; i<=N; i++) {
5077           F = Fmax*i/(double)N;
5078           k = (*model->k)(F, &env, params);
5079           if (k < 0) break;
5080           printf("%g\t%g\n", F, k);
5081         }
5082       }
5083       break;
5084     case M_SPECIAL :
5085       if (model->k == NULL || model->k == &bell_k) {
5086         printf("No special mode for model %s\n", model->name);
5087         exit(1);
5088       }
5089       if (special_xmax <= special_xmin) {
5090         printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5091         exit(1);
5092       }
5093       {
5094         double Xrng=(special_xmax-special_xmin), x, E;
5095         int i, N=500;
5096         printf("#x (m)\tE (J)\n");
5097         for (i=0; i<=N; i++) {
5098           x = special_xmin + Xrng*i/(double)N;
5099           E = multi_model_E(model->k, params, &env, x);
5100           printf("%g\t%g\n", x, E);
5101         }
5102       }
5103       break;
5104     default :
5105       break;
5106   }
5107   if (model->destructor)
5108     (*model->destructor)(params);
5109   return 0;
5110 }
5111 @
5112
5113 <<k model utility multi model E>>=
5114 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5115 {
5116   kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5117   double E;
5118   if (k == kramers_integ_k)
5119     E = (*p->E)(x, p->E_params);
5120   else
5121     E = kramers_E(0, env, model_params, x);
5122   return E;
5123 }
5124
5125 @
5126
5127 <<k model utility includes>>=
5128 #include <stdlib.h>
5129 #include <stdio.h>
5130 #include <string.h> /* strlen() */
5131 #include <assert.h> /* assert() */
5132 #include "global.h"
5133 #include "parse.h"
5134 #include "string_eval.h"
5135 #include "k_model.h"
5136 @
5137
5138 <<k model utility definitions>>=
5139 <<version definition>>
5140 #define VFLAG 1 /* verbose */
5141 enum mode_t {M_K_OF_F, M_SPECIAL};
5142 <<string comparison definition and globals>>
5143 <<k model getopt definitions>>
5144 <<initialize model definition>>
5145 <<kramers k structure>>
5146 <<kramers integ k structure>>
5147
5148 @
5149
5150 <<k model utility getopt functions>>=
5151 <<index k model>>
5152 <<k model utility help>>
5153 <<k model utility get options>>
5154 @
5155
5156 <<k model utility help>>=
5157 void help(char *prog_name,
5158           environment_t *env,
5159           int n_k_models, k_model_getopt_t *k_models,
5160           int k_model, double Fmax, double special_xmin, double special_xmax)
5161 {
5162   int i, j;
5163   printf("usage: %s [options]\n", prog_name);
5164   printf("Version %s\n\n", VERSION);
5165   printf("A utility for understanding the available unfolding models\n\n");
5166   printf("Environmental options:\n");
5167   printf("-T T\tTemperature T (currently %g K)\n", env->T);
5168   printf("-C T\tYou can also set the temperature T in Celsius\n");
5169   printf("Model options:\n");
5170   printf("-k model\tTransition rate model (currently %s)\n",
5171          k_models[k_model].name);
5172   printf("-K args\tTransition rate model argument string (currently %s)\n",
5173          k_models[k_model].params);
5174   printf("There are two output modes.  In standard mode, k(F) is printed\n");
5175   printf("For example:\n");
5176   printf("  #Force (N)\tk (% pop. per s)\n");
5177   printf("  123.456\t7.89\n");
5178   printf("  ...\n");
5179   printf("In special mode, the output depends on the model.\n");
5180   printf("For models defining an energy function E(x), that function is printed\n");
5181   printf("  #Distance (m)\tE (J)\n");
5182   printf("  0.0012\t0.0034\n");
5183   printf("  ...\n");
5184   printf("-m\tChange output to standard mode\n");
5185   printf("-M\tChange output to special mode\n");
5186   printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5187   printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5188   printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5189   printf("-V\tChange output to verbose mode\n");
5190   printf("-h\tPrint this help and exit\n");
5191   printf("\n");
5192   printf("Unfolding rate models:\n");
5193   for (i=0; i<n_k_models; i++) {
5194     printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5195     for (j=0; j < k_models[i].num_params ; j++)
5196       printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5197     printf("\t\tdefaults: %s\n", k_models[i].params);
5198   }
5199 }
5200 @
5201
5202 <<k model utility get options>>=
5203 void get_options(int argc, char **argv, environment_t *env,
5204                  int n_k_models, k_model_getopt_t *k_models,
5205                  enum mode_t *mode, k_model_getopt_t **model,
5206                  double *Fmax, double *special_xmin, double *special_xmax,
5207                  unsigned int *flags)
5208 {
5209   char *prog_name = NULL;
5210   char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5211   int k_model=0;
5212   extern char *optarg;
5213   extern int optind, optopt, opterr;
5214
5215   assert (argc > 0);
5216
5217   /* setup defaults */
5218
5219   prog_name = argv[0];
5220   env->T = 300.0;   /* K           */
5221   *mode = M_K_OF_F;
5222   *flags = 0;
5223   *model = k_models;
5224   *Fmax = 1e-9;
5225   *special_xmax = 1e-8;
5226   *special_xmin = 0.0;
5227
5228   while ((c=getopt(argc, argv, options)) != -1) {
5229     switch(c) {
5230     case 'T':  env->T   = atof(optarg);           break;
5231     case 'C':  env->T   = atof(optarg)+273.15;    break;
5232     case 'k':
5233       k_model = index_k_model(n_k_models, k_models, optarg);
5234       *model = k_models+k_model;
5235       break;
5236     case 'K':
5237       k_models[k_model].params = optarg;
5238       break;
5239     case 'm': *mode = M_K_OF_F;             break;
5240     case 'M': *mode = M_SPECIAL;            break;
5241     case 'F': *Fmax = atof(optarg);         break;
5242     case 'x': *special_xmin = atof(optarg); break;
5243     case 'X': *special_xmax = atof(optarg); break;
5244     case 'V': *flags |= VFLAG;              break;
5245     case '?':
5246       fprintf(stderr, "unrecognized option '%c'\n", optopt);
5247       /* fall through to default case */
5248     default:
5249       help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5250       exit(1);
5251     }
5252   }
5253   return;
5254 }
5255 @
5256
5257
5258 \section{\LaTeX\ documentation}
5259
5260 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5261 The comment blocks are extracted (with nicely formatted code blocks), using
5262 <<latex makefile lines>>=
5263 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5264         noweave -latex -index -delay $< > $@
5265 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5266         cp -f $< $@
5267 @
5268 and compiled using
5269 <<latex makefile lines>>=
5270 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5271                 | $(DOC_DIR)
5272         $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5273         $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5274         $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5275         $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5276         mv $(BUILD_DIR)/sawsim.pdf $@
5277 @
5278 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5279 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5280 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5281
5282 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5283 <<latex makefile lines>>=
5284 semi-clean_latex :
5285         rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5286                 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5287                 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
5288                 $(BUILD_DIR)/sawsim.bib
5289 clean_latex : semi-clean_latex
5290         rm -f $(DOC_DIR)/sawsim.pdf
5291 @
5292
5293 \section{Makefile}
5294
5295 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5296 In this case, we have several \emph{targets} that we'd like to build:
5297  the main simulation program \prog;
5298  the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5299  the documentation [[sawsim.pdf]];
5300  and the [[Makefile]] itself.
5301 Besides the generated files, there is the utility target
5302  [[clean]] for removing all generated files except [[Makefile]].
5303 % and
5304 % [[dist]] for generating a distributable tar file.
5305
5306 Extract the makefile with
5307 `[[notangle -Rmakefile src/sawsim.nw | sed 's/        /\t/' > Makefile]]'.
5308 The sed is needed because notangle replaces tabs with eight spaces,
5309 but [[make]] needs tabs.
5310 <<makefile>>=
5311 # Decide what directories to put things in
5312 SRC_DIR = src
5313 BUILD_DIR = build
5314 BIN_DIR = bin
5315 DOC_DIR = doc
5316 # Note: Cannot use, for example, `./build', because make eats the `./'
5317 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5318
5319 # Modules (X.c, X.h) defined in the noweb file
5320 NW_SAWSIM_MODS =
5321 CHECK_BINS =
5322 # Modules defined outside the noweb file
5323 FREE_SAWSIM_MODS = interp tavl
5324
5325 <<list module makefile lines>>
5326 <<tension balance module makefile lines>>
5327 <<tension model module makefile lines>>
5328 <<k model module makefile lines>>
5329 <<parse module makefile lines>>
5330 <<string eval module makefile lines>>
5331 <<accel k module makefile lines>>
5332
5333 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5334
5335 # Everything needed for sawsim under one roof.  sawsim.c must come first
5336 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5337         $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5338 # Libraries needed to compile sawsim
5339 LIBS = gsl gslcblas m
5340 CHECK_LIBS = gsl gslcblas m check
5341 # The non-check binaries generated
5342 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5343 DOCS = sawsim.pdf
5344
5345 # Define the major targets
5346 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5347
5348 view : $(DOC_DIR)/sawsim.pdf
5349         xpdf $< &
5350 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
5351         $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5352                 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5353         gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5354 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5355         $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5356 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5357                 clean_tension_model_utils \
5358                 clean_k_model_utils clean_latex clean_check_sawsim
5359         rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5360                 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5361                 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5362                 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5363                 $(BUILD_DIR)/global.h ./gmon.out
5364         $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5365
5366 # Various builds of sawsim
5367 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
5368         gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5369 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
5370         gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5371 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
5372         gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5373
5374 # Create the directories
5375 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5376         mkdir $@
5377
5378 # Copy over the source external to sawsim
5379 # Note: Cannot use, for example, `./build', because make eats the `./'
5380 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5381 .SECONDEXPANSION:
5382 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5383                 | $(BUILD_DIR)
5384         cp -f $< $@
5385 .SECONDEXPANSION:
5386 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5387                 | $(BUILD_DIR)
5388         cp -f $< $@
5389
5390 ## Basic source generated with noweb
5391 # The central files sawsim.c and global.h...
5392 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5393         notangle $< > $@
5394 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5395         notangle -Rglobal.h $< > $@
5396 # ... and the modules
5397 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5398         notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5399 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5400         notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5401 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5402         notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5403 # Note: I use `_' as a space in my file names, but noweb doesn't like
5404 # that so we use `-' to name the noweb roots and substitute here.
5405
5406 ## Building the unit-test programs
5407 # for sawsim.c...
5408 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5409         notangle -Rchecks $< > $@
5410 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5411                 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5412                 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
5413         gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5414 clean_check_sawsim :
5415         rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5416 # ... and the modules
5417 .SECONDEXPANSION:
5418 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5419                 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5420                 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5421                 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5422                 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5423                 $$(BUILD_DIR)/global.h | $(BIN_DIR)
5424         gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5425                 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5426                 $(CHECK_LIBS:%=-l%)
5427 # todo: clean up the required modules to
5428 $(CHECK_BINS:%=clean_%) :
5429         rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5430
5431 # Cleaning up the modules
5432 .SECONDEXPANSION:
5433 $(SAWSIM_MODS:%=clean_%) :
5434         rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5435
5436 <<tension model utils makefile lines>>
5437 <<k model utils makefile lines>>
5438 <<latex makefile lines>>
5439
5440 Makefile : $(SRC_DIR)/sawsim.nw
5441         notangle -Rmakefile $< | sed 's/        /\t/' > $@
5442 @
5443 This is fairly self-explanatory.  We compile a dynamically linked
5444 version ([[sawsim]]) and a statically linked version
5445 ([[sawsim_static]]).  The static version is about 9 times as big, but
5446 it runs without needing dynamic GSL libraries to link to, so it's a
5447 better format for distributing to a cluster for parallel evaluation.
5448
5449 \section{Math}
5450
5451 \subsection{Hookean springs in series}
5452 \label{app.math_hooke}
5453
5454 \begin{theorem}
5455 The effective spring constant for $n$ springs $k_i$ in series is given by
5456 $$
5457   k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5458 $$
5459 \end{theorem}
5460
5461 \begin{proof}
5462 \begin{align}
5463   F &= k_i x_i &&\forall i \le n \\
5464   x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5465   x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5466   F &= k_1 x_1 = k_\text{eff} x \\
5467   k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5468                 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5469 \end{align}
5470 \end{proof}
5471
5472 \phantomsection
5473 \addcontentsline{toc}{section}{References}
5474 \bibliography{sawsim}
5475
5476 \end{document}
5477