Fixed 'y > yub' error by reducing default tolerances min_d{x,y} in

[sawsim.git] / sawsim.nw

%%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
% we have our own preamble, 
% use `noweave -delay` to not print noweb's automatic one
  % -*- mode: Noweb; noweb-code-mode: c-mode -*-
\documentclass[letterpaper, 11pt]{article}
\usepackage{noweb}
\pagestyle{noweb}
\noweboptions{smallcode}

\usepackage{url}    % needed for natbib for underscores in urls?
\usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
% breaklinks breaks long links
% colorlinks colors link text instead of boxing it
\usepackage{amsmath}  % for \text{}, \xrightarrow[sub]{super}
\usepackage{amsthm}  % for theorem and proof environments
\newtheorem{theorem}{Theorem}

\usepackage{doc}    % BibTeX
\usepackage[super,sort&compress]{natbib} % fancy citation extensions
% super selects citations in superscript mode
% sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})

%\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
%\bibliographystyle{plain} % pick the bibliography style, includes dates
\bibliographystyle{plainnat}
\defcitealias{sw:noweb}{{\tt noweb}}
\defcitealias{galassi05}{GSL}
\defcitealias{sw:check}{{\tt check}}
\defcitealias{sw:python}{Python}

\topmargin -1.0in
\headheight 0.5in
\headsep 0.5in
\textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
\oddsidemargin -0.5in
\textwidth 7.5in
\setlength{\parindent}{0pt}
\setlength{\parskip}{5pt}

% For some reason, the \maketitle wants to go on it's own page
% so we'll just hardcode the following in our first page.
%\title{Sawsim: a sawtooth protein unfolding simulator}
%\author{W.~Trevor King}
%\date{\today}

\newcommand{\prog}{[[sawsim]]}
\newcommand{\lang}{[[C]]}
\newcommand{\p}[3]{\left#1 #2 \right#3}   % e.g. \p({complicated expr})
\newcommand{\Th}{\ensuremath{\sup{\text{th}}}}   % ordinal th
\newcommand{\U}[1]{\ensuremath{\text{ #1}}}      % units: $1\U{m/s}$

% single spaced lists, from Alvin Alexander
% http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
\newenvironment{packed_item}{
\begin{itemize}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{itemize}}

\begin{document}
\nwfilename{sawsim.nw}
\nwbegindocs{0}

%\maketitle
\begin{centering}
  Sawsim: a sawtooth protein unfolding simulator \\
  W.~Trevor King \\
  \today \\
\end{centering}
\bigskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

The unfolding of multi-globular protein strings is a stochastic
process.  In the \prog\ program, we use Monte Carlo methods to
simulate this unfolding through an explicit time-stepping approach.
We develop a program in \lang\ to simulate probability of unfolding
under a constant extension velocity of a chain of globular domains.

\subsection{Related work}

TODO References

\subsection{About this document}

This document was written in \citetalias{sw:noweb} with code blocks in
\lang\ and comment blocks in \LaTeX.

\subsection{Change Log}

Version 0.0 used only the unfolded domain WLC to determine the tension
and had a fixed timestep $dt$.

Version 0.1 added adaptive timesteps, adjusting so that the unfolding
probability for a given domain was constant.

Version 0.2 added Kramers' model unfolding rate calculations, and a
switch to select Bell's or Kramers' model.  This induced a major
rewrite, introducing the tension group abstraction, and a split to
multiple \lang\ source files.  Also added a unit-test suites for
sanity checking, and switched to SI units throughout.

Version 0.3 added integrated Kramers' model unfolding rate
calculations.  Also replaced some of my hand-coded routines with GNU
Scientific Library \citepalias{galassi05} calls.

Version 0.4 added Python evaluation for the saddle-point Kramers'
model unfolding rate calculations, so that the model functions could
be controlled from the command line.  Also added interpolating lookup
tables to accelerate slow unfolding rate model evaluation, and command
line control of tension groups.

<<version definition>>=
#define VERSION "0.4"
@ 

\subsection{License}

<<license>>=
 sawsim - program for simulating single molecule mechanical unfolding.
 Copyright (C) 2008, William Trevor King

 This program is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public License as
 published by the Free Software Foundation; either version 3 of the
 License, or (at your option) any later version.

 This program is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
 See the GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 02111-1307, USA.

 The author may be contacted at <wking@drexel.edu> on the Internet, or
 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
 Philadelphia PA 19104, USA.
@ 

<<license comment>>=
/*
 <<license>>
 */

@ 

\section{Main}

The unfolding system is stored as a chain of \emph{domains}.  Domains
can be folded, globular protein domains, unfolded protein linkers, AFM
cantilevers, or other stretchable system components.

Each domain contributes to the total tension, which depends on the
chain extension.  The particular model for the tension as a function
of extension is handled generally with function pointers.  So far the
following models have been implemented:
\begin{packed_item}
 \item  Null (Appendix \ref{sec.null_tension}),
 \item  Constant (Appendix \ref{sec.const_tension}),
 \item  Hookean spring (Appendix \ref{sec.hooke}),
 \item  Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
 \item  Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
\end{packed_item}

The tension model and parameters used for a particular domain depend
on whether the domain is folded or unfolded.  The transition rate from
the folded state to the unfolded state is also handled generally with
function pointers, with implementations for
\begin{packed_item}
 \item  Null (Appendix \ref{sec.null_k}),
 \item  Bell's model (Appendix \ref{sec.bell}),
 \item  Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
 \item  Kramers' saddle point model (Appendix \ref{sec.kramers}).
\end{packed_item}

The unfolding of the chain domains is modeled in two stages.  First
the tension of the chain is calculated.  Then the tension (assumed to
be constant over the length of the chain), is applied to each folded
domain, and used to calculate the probability of that subdomain
unfolding in the next timestep $dt$.  Then the time is stepped
forward, and the process is repeated.

<<main program>>=
int main(int argc, char **argv)
{
  <<tension model getopt array definition>>
  <<k model getopt array definition>>
  list_t *domain_list=NULL;    /* variables initialized in get_options() */
  environment_t env={0};
  double P, dt_max, v;
  int num_folded=0, unfolding;
  double x=0, dt, F;           /* variables used in the simulation loop */
  get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_MODELS,
            tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
  setup();
  if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n");
  if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
  while (num_folded > 0) {
    F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
    dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
    unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
    if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
    x += v*dt;
  }
  destroy_domain_list(domain_list);
  free_accels();
  return 0;
}
@ The meat of the simulation is bundled into the three functions
[[find_tension]], [[determine_dt]], and [[random_unfoldings]].
[[find_tension]] is discussed in Section \ref{sec.find_tension},
[[determine_dt]] in Section \ref{sec.adaptive_dt}, and
[[random_unfoldings]] in Section \ref{sec.unfolding_rate}.

Environmental parameters important in determining reaction rates and
tensions (e.g. temperature, pH) are stored in a single structure to
facilitate extension to more complicated models in the future.
<<environment definition>>= 
typedef struct environment_struct {
  double T;
} environment_t;
@ 

<<global.h>>=
<<environment definition>>
<<one dimensional function definition>>
<<create/destroy definitions>>
<<model globals>>
@ 

\section{Simulation functions}

<<simulation functions>>=
<<find tension>>
<<determine dt>>
<<happens>>
<<does domain unfold>>
<<randomly unfold domains>>
@ 

\subsection{Tension}
\label{sec.find_tension}

Because the stretched system may be made up of several parts (folded
domains, unfolded domains, spring-like cantilever, \ldots), we will
assign the domains to models and groups.  The models are used to
determine a domain's tension (Hookean spring, WLC, \ldots).  Groups
will represent collections of domains which the model should treat as
a single entity.  A domain may belong to different groups or models
depending on its state.  For example, a domain might be modeled as a
freely-jointed chain when it iss folded and as a worm-like chain class
when it is unfolded.  The domains are assumed to be commutative, so
ordering is ignored.  The interactions between the groups are assumed
to be linear, but the interactions between domains of the same group
need not be.  This allows for non-linear group models such as th
worm-like or freely-jointed chains.  Each model has a tension handler
function, which gives the tension $F$ needed to stretch a given group
of domains a distance $x$.

To understand the purpose of groups, consider a chain of proteins
where two folded proteins $A$ and $B$ are modeled as WLCs with
persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
modeled as WLCs with persistence length $p_u$.  The proteins are
attached to a cantilever $E$ qof spring constant $κ$.  The string
would get split into three lists:
\begin{center}
  \begin{tabular}{llll}
    Model & Group &    List & Tension \\
    WLC   &     0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B}$ \\
    WLC   &     1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D}$ \\
    Hooke &     0 & $[E]$   & $F_\text{Hooke}(x, κ}$ \\
  \end{tabular}
\end{center}
Note that group numbers only matter \emph{within} model classes, since
grouping domains with seperate models doesn't make sense.

<<find tension definitions>>=
#define NUM_TENSION_MODELS 5

@ 

<<tension handler array global declaration>>=
extern one_dim_fn_t *tension_handlers[];
@

<<tension handler array global>>= 
one_dim_fn_t *tension_handlers[] = {
              NULL,
              &const_tension_handler,
              &hooke_handler,
              &wlc_handler,
              &fjc_handler,
              };

@ 

<<tension model global declarations>>=
<<tension handler array global declaration>>
@ 

<<tension handler types>>=
typedef struct tension_handler_data_struct {
  /* int            num_domains; */
  /* domain_t      *domains;     */
  list_t        *group;
  environment_t *env;
  void          *persist;
} tension_handler_data_t;

@ 

After sorting the chain into separate groups $G_i$, with tension
handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
\forall i,j$.  Note that there may be several groups within each model
class.  [[find_tension]] is basically a wrapper to feed proper input
into the [[tension_balance]] function.  [[unfolding]] should be set to
0 if there were no unfoldings since the previous call to
[[find_tension]].
<<find tension>>=
double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
{
  static int num_active_groups=0;
  static one_dim_fn_t **per_group_handlers = NULL;
  static void **per_group_params = NULL;
  static double last_x;
  tension_handler_data_t *tdata;
  double F;
  int i;

#ifdef DEBUG
  fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
  list_print(stderr, domain_list, "domain list");
#endif

  if (unfolding != 0 || num_active_groups == 0) { /* setup statics */  
    /* free old statics */
    if (per_group_handlers != NULL)
      free(per_group_handlers);
    if (per_group_params != NULL) {
      for (i=0; i < num_active_groups; i++) {
        assert(per_group_params[i] != NULL);
        free(per_group_params[i]);
      }
      free(per_group_params);
    }

    /* get new statics */
    get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);

    /* fill in the group handlers and params */
    for (i=0; i<num_active_groups; i++) {
      tdata = (tension_handler_data_t *) per_group_params[i];
#ifdef DEBUG
      fprintf(stderr, "active group %d of %d", i, num_active_groups);
      list_print(stderr, tdata->group, " parameters");
#endif
      tdata->env = env;
      /* tdata->persist continues unchanged */
    }
    last_x = -1.0;
  } /* else roll with the current statics */

  F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
  last_x = x;
  return F;
}
@ For the moment, we will restrict our group boundaries to a single
dimension, so $\sum_i x_i = x$, our total extension (it is this
restriction that causes the groups to interact linearly).  We'll also
restrict our extensions to all be positive.  With these restrictions,
the problem of balancing the tensions reduces to solving a system of
$N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
the number of active groups.  In general this can be fairly
complicated, but without much loss of practicality we can restrict
ourselves to strictly monotonically increasing, non-negative tension
functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
simpler.  For example, it guarantees the existence of a unique, real
solution for finite forces.  See Appendix \ref{app.tension_balance}
for the tension balancing implementation.

Each group must define a parameter structure if the tension model is
parametrized,
 a function to create the parameter structure,
 a function to destroy the parameter structure, and
 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
The tension-specific parameter structures are contained in the domain
group field.  See the Section \ref{app.model_selection} for a
discussion on the command line framework.  See the worm-like chain
implementation for an example (Section \ref{sec.wlc}).

The major limitation of our approach is that all of our tension models
are Markovian (memory-less), which is not really a problem for the
chains (see \citet{evans99} for WLC thermalization rate calculations).

\subsection{Unfolding rate}
\label{sec.unfolding_rate}

Each folded domain is modeled as unimolecular, first order reaction
$$
  \text{Folded} \xrightarrow{k}  % in tex: X atop Y
  \text{Unfolded}
$$
With the rate constant $k$ defined by
$$
  \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
$$
So the probability of a given protein unfolding in a short time $dt$
is given by
$$
  P = k \cdot dt.
$$

<<does domain unfold>>=
int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
  double k;
  k = accel_k(domain->k, F, env, domain->k_params);
  //(*domain->k)(F, env, domain->k_params);
  //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
  return happens(k*dt); /* dice roll for prob. k*dt event */
}
@ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.

Only one domain can unfold in each timestep, because the timescale of
a domain unfolding $dt_u$ is assumed to be less than the simulation
timestep $dt$, so a domain will completely unfold in a single
timestep.  We adapt our timesteps to keep the probability of a single
domain unfolding low, so the probability of two domains unfolding in
the same timestep is negligible.  This approach breaks down as the
adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
shouldn't be a problem.  To reassure yourself, you can ask the
simulator to print the smallest timestep that was used with TODO.
<<randomly unfold domains>>=
int random_unfoldings(list_t *domain_list, double tension, 
                      double dt, environment_t *env,
                      int *num_folded_domains)
{ /* returns 1 if there was an unfolding and 0 otherwise */
  while (domain_list != NULL) {
    if (D(domain_list)->state == DS_FOLDED 
        && domain_unfolds(tension, dt, env, D(domain_list))) {
      if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
        fprintf(stdout, "%g\n", tension);
      D(domain_list)->state = DS_UNFOLDED;
      (*num_folded_domains)--;
      return 1; /* our one domain has unfolded, stop looking for others */
    }
    domain_list = domain_list->next;
  }
  return 0;
}
@ 

\subsection{Adaptive timesteps}
\label{sec.adaptive_dt}

We'd like to pick $dt$ so the probability of unfolding in the next
timestep is small.  If we don't adapt our timestep, we risk breaking
our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
only one domain unfolds in a given timestep.  Because $F(x)$ is
monotonically increasing, excessively large timesteps will lead to
erroneously large unfolding forces.  The simulation would have been
accurate for sufficiently small fixed timesteps, but adaptive
timesteps will allow us to move through low tension regions in fewer
steps, leading to a more efficient simulation.

The actual adaptive timestep implementation is not particularly
interesting, since we are only required to reduce $dt$ to somewhere
below a set threshold, so I've removed it to Appendix
\ref{app.adaptive_dt}.

\section{Models}
 
TODO: model intro.

The models provide the physics of the simulation, but the simulation
allows interchangeable models, and we are currently focusing on the
simulation itself, so we remove the actual model implementation to the
appendices.

The tension models are defined in Appendix \ref{sec.tension_models}
and unfolding models are defined in Appendix \ref{sec.k_models}.

<<model globals>>=
#define k_B   1.3806503e-23  /* J/K */
@ 


\section{Command line interface}

<<get option functions>>=
<<help>>
<<index tension model>>
<<index k model>>
<<generate domain>>
<<get options>>
@ 

\subsection{Model selection}
\label{app.model_selection}

The main difficulty with the command line interface in \prog\ is
developing an intuitive method for accessing the possibly large number
of available models.  We'll treat this generally by defining an array
of available models, containing their important parameters.

<<get options definitions>>=
<<model getopt definitions>>
@ 

<<create/destroy definitions>>=
typedef void *create_data_func_t(char **param_strings);
typedef void destroy_data_func_t(void *param_struct);
@ 

<<model getopt definitions>>=
<<tension model getopt definitions>>
<<k model getopt definitions>>
@ 


\subsubsection{Tension}

To access the various tension models from the command line, we define
a structure that contains all the neccessary information about the
model.
<<tension model getopt definitions>>=
typedef struct tension_model_getopt_struct {
  one_dim_fn_t      *handler;     /* fn implementing the model on a group */
  char                *name;        /* name string identifying the model    */
  char                *description; /* a brief description of the model     */
  int                  num_params;  /* number of extra params for the model */
  char               **param_descriptions; /* descriptions of extra params  */
  char                *params;      /* default values for the extra params  */
  create_data_func_t  *creator;     /* param-string -> model-data-struct fn */
  destroy_data_func_t *destructor;  /* fn to free the model data structure  */
} tension_model_getopt_t;
@ 

<<tension model getopt array definition>>=
tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
  <<null tension model getopt>>,
  <<constant tension model getopt>>,
  <<hooke tension model getopt>>,
  <<worm-like chain tension model getopt>>,
  <<freely-jointed chain tension model getopt>>
};
@ 

\subsubsection{Unfolding rate}

<<k model getopt definitions>>=
#define NUM_K_MODELS 5

typedef struct k_model_getopt_struct {
  char *name;
  char *description;
  k_func_t *k;
  int num_params;
  char **param_descriptions;
  char *params;
  create_data_func_t *creator;
  destroy_data_func_t *destructor;
} k_model_getopt_t;
@ 

<<k model getopt array definition>>=
k_model_getopt_t k_models[NUM_K_MODELS] = {
  <<null k model getopt>>,
  <<const k model getopt>>,
  <<bell k model getopt>>,
  <<kramers k model getopt>>,
  <<kramers integ k model getopt>>
};
@ 

\subsection{help}

<<help>>=
void help(char *prog_name, double P, double dt_max, double v,
          environment_t *env,
          int n_tension_models, tension_model_getopt_t *tension_models,
          int folded_tension_model, int unfolded_tension_model,
          int n_k_models, k_model_getopt_t *k_models,
          int k_model)
{
  int i, j;
  printf("usage: %s [options]\n", prog_name);
  printf("Version %s\n\n", VERSION);
  printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
  printf("Simulation options:\n");
  printf("-P P\tTarget probability for dt (currently %g)\n", P);
  printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
  printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
  printf("Environmental options:\n");
  printf("-T T\tTemperature T (currently %g K)\n", env->T);
  printf("-C T\tYou can also set the temperature T in Celsius\n");
  printf("Model options:\n");
  printf("The domains exist in either the folded or unfolded state\n");
  printf("The following options change the default behavior in each state.\n");
  printf("-m model[,group]\tFolded tension model (currently %s)\n",
         tension_models[folded_tension_model].name);
  printf("-a args\tFolded tension model argument string (currently %s)\n",
         tension_models[folded_tension_model].params);
  printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
         tension_models[unfolded_tension_model].name);
  printf("-A args\tUnfolded tension model argument string (currently %s)\n",
         tension_models[unfolded_tension_model].params);
  printf("The following options change the unfolding rate.\n");
  printf("-k model\tTransition rate model (currently %s)\n",
         k_models[k_model].name);
  printf("-K args\tTransition rate model argument string (currently %s)\n",
         k_models[k_model].params);
  printf("Domain creation options:\n");
  printf("Once you've set up the appropriate default models, you need to add the domains\n");
  printf("-F n\tAdd n folded domains with the current models\n");
  printf("-U n\tAdd n unfolded domains with the current models\n");
  printf("Output mode options:\n");
  printf("There are two output modes.  In standard mode, only the unfolding\n");
  printf("events are printed.  For example:\n");
  printf("  #Unfolding Force (N)\n");
  printf("  123.456e-12\n");
  printf("  ...\n");
  printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
  printf("  #Position (m)\tForce (N)\n");
  printf("  0.001\t0.0023\n");
  printf("  ...\n");
  printf("-V\tChange output to verbose mode\n");
  printf("-h\tPrint this help and exit\n");
  printf("\n");
  printf("Tension models:\n");
  for (i=0; i<n_tension_models; i++) {
    printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
    for (j=0; j < tension_models[i].num_params ; j++)
      printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
    printf("\t\tdefaults: %s\n", tension_models[i].params);
  }
  printf("Unfolding rate models:\n");
  for (i=0; i<n_k_models; i++) {
    printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
    for (j=0; j < k_models[i].num_params ; j++)
      printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
    printf("\t\tdefaults: %s\n", k_models[i].params);
  }
  printf("Examples:\n");
  printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
  printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
  printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
  printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
}
@ 

\subsection{Get options}

<<get options>>=
void get_options(int argc, char **argv,
                 double *pP, double *pDt_max, double *pV, 
                 environment_t *env,
                 int n_tension_models, tension_model_getopt_t *tension_models,
                 int n_k_models, k_model_getopt_t *k_models,
                 list_t **pDomain_list, int *pNum_folded)
{
  char *prog_name = NULL;
  char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh";
  int ftension_model=0, utension_model=0, k_model=0;
  char *ftension_params=NULL, *utension_params=NULL;
  int i, n, ftension_group=0, utension_group=0;
  int num_strings;
  char **strings;
  extern char *optarg;
  extern int optind, optopt, opterr;

  assert (argc > 0);

#ifdef DEBUG
  for (i=0; i<n_tension_models; i++) {
    fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
            i, tension_models[i].name, tension_models[i].handler);
    assert(tension_models[i].handler == tension_handlers[i]);
  }
#endif

  /* setup defaults */
  flags = FLAG_OUTPUT_UNFOLDING_FORCES;
  prog_name = argv[0];
  *pP = 1e-3;       /* % pop per s */
  *pDt_max = 0.001; /* s           */
  *pV = 1e-6;       /* m/s         */
  env->T = 300.0;   /* K           */

  while ((c=getopt(argc, argv, options)) != -1) {
    switch(c) {
    case 'P':  *pP      = atof(optarg);           break;
    case 't':  *pDt_max = atof(optarg);           break;
    case 'v':  *pV      = atof(optarg);           break;
    case 'T':  env->T   = atof(optarg);           break;
    case 'C':  env->T   = atof(optarg)+273.15;    break;
    case 'm':
      parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
      assert(num_strings == 1 || num_strings == 2);
      if (num_strings == 1)
        ftension_group = 0;
      else
        ftension_group = atoi(strings[1]);
      ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
      ftension_params = tension_models[ftension_model].params;
      free_parsed_list(num_strings, strings);
      break;
    case 'a':
      ftension_params = optarg;
      break;
    case 'M':
      parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
      assert(num_strings == 1 || num_strings == 2);
      if (num_strings == 1)
        utension_group = 0;
      else
        utension_group = atoi(strings[1]);
      utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
      utension_params = tension_models[utension_model].params;
      free_parsed_list(num_strings, strings);
      break;
    case 'A':
      utension_params = optarg;
      break;
    case 'k':
      k_model = index_k_model(n_k_models, k_models, optarg);
      break;
    case 'K':
      k_models[k_model].params = optarg;
      break;
    case 'F':
      n = atoi(optarg);
      assert(n > 0);
      for (i=0; i<n; i++) {
        push(pDomain_list, generate_domain(DS_FOLDED,
                                           tension_models+ftension_model,
                                           ftension_group,
                                           ftension_params,
                                           tension_models+utension_model,
                                           utension_group,
                                           utension_params,
                                           k_models+k_model), 1);
      }
      *pNum_folded += n;
      break;
    case 'U':
      n = atoi(optarg);
      assert(n > 0);
      for (i=0; i<n; i++) {
        push(pDomain_list, generate_domain(DS_UNFOLDED,
                                           tension_models+ftension_model,
                                           ftension_group,
                                           ftension_params,
                                           tension_models+utension_model,
                                           utension_group,
                                           utension_params,
                                           k_models+k_model), 1);
      }
      break;
    case 'V':
      flags = FLAG_OUTPUT_FULL_CURVE;
      break;
    case '?':
      fprintf(stderr, "unrecognized option '%c'\n", optopt);
      /* fall through to default case */
    default:
      help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
      exit(1);
    }
  }
  /* check the arguments */
  assert(*pDomain_list != NULL); /* no domains! */
  assert(*pP > 0.0); assert(*pP < 1.0);
  assert(*pDt_max > 0.0);
  assert(*pV > 0.0);
  assert(env->T > 0.0);
  return;
}
@ 

<<index tension model>>=
int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
{
  int i;
  for (i=0; i<n_models; i++)
    if (STRMATCH(models[i].name, name))
      break;
  if (i == n_models) {
    fprintf(stderr, "Unrecognized tension model '%s'\n", name);
    exit(1);
  }
  return i;
}
@ 

<<index k model>>=
int index_k_model(int n_models, k_model_getopt_t *models, char *name)
{
  int i;
  for (i=0; i<n_models; i++)
    if (STRMATCH(models[i].name, name))
      break;
  if (i == n_models) {
    fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
    exit(1);
  }
  return i;
}
@ 

<<initialize model definition>>=
/* requires int num_param_args and char **param_args in the current scope
 * usage:
 *  INIT_MODEL("folded", folded_model, folded_params);
 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
 */
#define INIT_MODEL(role, model, param_string, param_pointer) \
  do {                                                       \
    parse_list_string(param_string, SEP, '{', '}',           \
                      &num_param_args, &param_args);         \
    if (num_param_args != model->num_params) {               \
      fprintf(stderr,                                        \
              "%s model %s expected %d params,"              \
              role, model->name, model->num_params);         \
      fprintf(stderr,                                        \
              "not the %d params in '%s'\n",                 \
              num_param_args, param_string);                 \
      assert(num_param_args == model->num_params);           \
    }                                                        \
    if (model->creator)                                      \
      param_pointer = (*model->creator)(param_args);         \
    else                                                     \
      param_pointer = NULL;                                  \
    free_parsed_list(num_param_args, param_args);            \
  } while (0);
@ 

<<generate domain>>=
void *generate_domain(enum domain_state_t state,
                      tension_model_getopt_t *folded_model,
                      int folded_group,
                      char *folded_param_string,
                      tension_model_getopt_t *unfolded_model,
                      int unfolded_group,
                      char *unfolded_param_string,
                      k_model_getopt_t *k_model)
{
  void *folded_params, *unfolded_params, *k_params;
  int num_param_args; /* for INIT_MODEL() */
  char **param_args;  /* for INIT_MODEL() */

#ifdef DEBUG
  fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
  fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
          k_model->name, k_model->params,
          folded_model->name, folded_model->handler, folded_group, folded_param_string,
          unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
#endif
  INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
  INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
  INIT_MODEL("k", k_model, k_model->params, k_params);
  return create_domain(state,
                       k_model->k, k_params, k_model->destructor,
                       folded_model->handler, folded_group, folded_params, folded_model->destructor,
                       unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
                       );
}
@ 

\phantomsection
\appendix
\addcontentsline{toc}{section}{Appendicies}

\section{sawsim.c details}

\subsection{Layout}

The general layout of our simulation code is:
<<*>>=
<<license comment>>
<<includes>>
<<definitions>>
<<globals>>
<<functions>>
<<main program>>
@ 

We include [[math.h]], so don't forget to link to the libm with `-lm'.
<<includes>>=
#include <assert.h> /* assert()                                */
#include <stdlib.h> /* malloc(), free(), rand()                */
#include <stdio.h>  /* fprintf(), stdout                       */
#include <string.h> /* strlen, strtok()                        */
#include <math.h>   /* exp(), M_PI, sqrt()                     */
#include <sys/types.h> /* pid_t (returned by getpid())         */
#include <unistd.h> /* getpid() (for seeding rand()), getopt() */
#include "global.h"
#include "list.h"
#include "tension_balance.h"
#include "k_model.h"
#include "tension_model.h"
#include "parse.h"
#include "accel_k.h"
@ 

<<definitions>>=
<<version definition>>
<<flag definitions>>
<<max/min definitions>>
<<string comparison definition and globals>>
<<initialize model definition>>
<<get options definitions>>
<<domain definitions>>
@ 

<<globals>>=
<<flag globals>>
<<model globals>>
@ 

<<functions>>=
<<create/destroy domain>>
<<destroy domain list>>
<<group functions>>
<<simulation functions>>
<<boilerplate functions>>
@ 

<<boilerplate functions>>=
<<setup>>
<<get option functions>>
@ 

<<setup>>=
void setup(void)
{
  srand(getpid()*time(NULL)); /* seed rand() */
}
@ 

<<flag definitions>>=
/* in octal b/c of prefixed '0' */
#define FLAG_OUTPUT_FULL_CURVE       01
#define FLAG_OUTPUT_UNFOLDING_FORCES 02
@

<<flag globals>>=
static unsigned long int flags = 0;
@

\subsection{Utilities}
\label{app.utils}

<<max/min definitions>>=
#define MAX(a,b) ((a)>(b) ? (a) : (b))
#define MIN(a,b) ((a)<(b) ? (a) : (b))
@ 

Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
<<string comparison definition and globals>>=
// Check if two strings match, return 1 if they do
static char *temp_string_A;
static char *temp_string_B;
#define STRMATCH(a,b)   (temp_string_A=a, temp_string_B=b, \
        strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
        !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
        /* +1 to also compare the '\0' */
@ 

We also define a macro for our [[check]] unit testing
<<check relative error>>=
#define CHECK_ERR(max_err, expected, received)                               \
  do {                                                                       \
    fail_unless( (received-expected)/expected < max_err,                     \
                 "relative error %g >= %g in %s (Expected %g, received %g)", \
                 (received-expected)/expected, max_err, #received,           \
                 expected, received);                                        \
    fail_unless(-(received-expected)/expected < max_err,                     \
                 "relative error %g >= %g in %s (Expected %g, received %g)", \
                -(received-expected)/expected, max_err, #received,           \
                 expected, received);                                        \
  } while(0)
@ 

<<happens>>=
int happens(double probability)
{
  assert(probability >= 0.0); assert(probability <= 1.0);
  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*/
}
@


\subsection{Adaptive timesteps}
\label{app.adaptive_dt}

$F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
so basing the timestep on the the unfolding probability at the current tension
is dangerous, and we need to search for a $dt$ for which
$P(F(x+v*dt)) < P_\text{target}$.
There are two cases to consider.
In the most common, no domains have unfolded since the last step,
and we expect the next step to be slightly shorter than the previous one.
In the less common, domains did unfold in the last step,
and we expect the next step to be considerably longer than the previous one.
<<search dt>>=
double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
                 list_t *domain_list, 
                 environment_t *env, double x,
                 double target_prob, double max_dt, double v)
{ /* Returns the timestep dt in seconds for the current folded domain.
   * Takes a list of tension handlers, the list of domains,
   * a pointer env to the environmental data, a starting separation x in m,
   * a target_prob between 0 and 1,
   * max_dt in s, stretching velocity v in m/s.
   */
  double F, k, dtCur, dtU, dtUCur, dtL, dt;

  /* get upper bound using the current position */
  F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
  //printf("Start with x = %g (F = %g)\n", x, F);
  k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
  //printf("x %g\tF %g\tk %g\n", x, F, k);
  dtU = target_prob / k;    /* P = k dt, dtU is an upper bound on dt */
  if (dtU > max_dt) {
    //printf("overshot max_dt\n");  
    dtU = max_dt;
  }
  /* set a lower bound on dt too */
  dtL = 0.0;

  /* The dt determined above may produce illegitimate forces or ks.
   * Reduce the upper bound until we have valid ks. */
  dt = dtU;
  F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
  while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
    dtU /= 2.0;
    dt = dtU;
    F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
  }
  //printf("Try for dt = %g (F = %g)\n", dt, F);
  k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
  /* returned k may be -1.0 */
  //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); 
  while (k == -1.0) { /* reduce step until we hit a valid k */
    dtU /= 2.0;
    dt = dtU; /* hopefully, we can use the max dt, see if it works */
    F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
    //printf("Try for dt = %g (F = %g)\n", dt, F);
    k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
    //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); 
  }
  assert(dtU > 1e-14);      /* timestep to valid k too small */
  dtUCur = target_prob / k; /* safe timestep back from x+dtU */
  if (dtUCur >= dt)
    return dt; /* dtU is safe. */

  /* dtU wasn't safe, lets see what would be. */
  while (dtU > 1.1*dtL) { /* until the range is reasonably small */
    dt = (dtU + dtL) / 2.0;
    F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
    //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
    k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
    dtCur = target_prob / k;
    //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
    if (dtCur > dt) /* safe timestep back from x+dt covers dt */
      dtL = dt;
    else if (dtCur < dt) {  /* unsafe timestep back from x+dt, but...    */
      dtU = dt;             /* ... stepping out only dtCur would be safe */ 
      dtUCur = dtCur;
    } else
      break; /* dtCur = dt */
  }
  return MAX(dtUCur, dtL);
}
@ 

To determine $dt$ for an array of potentially different folded domains,
we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
<<determine dt>>=
<<search dt>>
double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
                    environment_t *env, double x,
                    double target_prob, double dt_max, double v)
{ /* Returns the timestep dt in seconds.
   * Takes the list of folded domains, target_prob between 0 and 1,
   * F in N, and T in K. */
  double dt=dt_max, new_dt;
  assert(target_prob > 0.0); assert(target_prob < 1.0);
  assert(dt_max > 0.0);
  
  /* .5 nm steps = v * dt */
  //return 0.5e-9/v;
  while (domain_list != NULL) {
    if (D(domain_list)->state == DS_FOLDED) {
      new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
      dt = MIN(dt, new_dt);
    }
    domain_list = domain_list->next;
  }
  return dt;
}
@ 

\subsection{Domain data}

Currently domains exist in two states, folded and unfolded, and the
only allowed transitions are folded $\rightarrow$ unfolded.  Of
course, it wouldn't be too complicated to extent this to a multi-state
system, with an array containing the domains group for each possible
state, and a matrix of transition-rate-calculating functions.
However, at this point such generality seems unnecessary at this
point.
<<domain definitions>>=
enum domain_state_t {DS_FOLDED,
                     DS_UNFOLDED
};

typedef struct domain_struct {
  enum domain_state_t state;
  one_dim_fn_t       *folded_handler;
  int                   folded_group;
  one_dim_fn_t       *unfolded_handler;
  int                   unfolded_group;
  k_func_t *k;           /* function returning unfolding rate */
  void *folded_params;   /* pointer to folded parameters   */
  void *unfolded_params; /* pointer to unfolded parameters */
  void *k_params;        /* pointer to k parameters        */
  destroy_data_func_t *destroy_folded;
  destroy_data_func_t *destroy_unfolded;
  destroy_data_func_t *destroy_k;
} domain_t;

/* get the domain data for the current list node */
#define D(list) ((domain_t *)(list)->d)
/* get the tension params for the current list node */
#define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
                     ? ((domain_t *)(list)->d)->folded_params   \
                     : ((domain_t *)(list)->d)->unfolded_params)
@ 
[[k]] is a pointer to the function determining the unfolding rate for a given tension.
[[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
[[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
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.

[[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
<<create/destroy domain>>=
domain_t *create_domain(enum domain_state_t state,
                        k_func_t *k,
                        void *k_params,
                        destroy_data_func_t *destroy_k,
                        one_dim_fn_t *folded_handler,
                        int folded_group,
                        void *folded_params,
                        destroy_data_func_t *destroy_folded,
                        one_dim_fn_t *unfolded_handler,
                        int unfolded_group,
                        void *unfolded_params,
                        destroy_data_func_t *destroy_unfolded)
{
  domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
  assert(ret != NULL);
  if (state == DS_FOLDED) {
    assert(k != NULL);  /* the pointer points somewhere valid  */
    assert(*k != NULL); /* and there is something useful there */
  } else
    assert(state == DS_UNFOLDED);
  ret->state = state;
  ret->folded_handler = folded_handler;
  ret->folded_group = folded_group;
  ret->unfolded_handler = unfolded_handler;
  ret->unfolded_group = unfolded_group;
  ret->k = k;
  ret->folded_params = folded_params;
  ret->unfolded_params = unfolded_params;
  ret->k_params = k_params;
  ret->destroy_folded = destroy_folded;
  ret->destroy_unfolded = destroy_unfolded;
  ret->destroy_k = destroy_k;
#ifdef DEBUG
  fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
#endif
  return ret;
}

void destroy_domain(domain_t *domain)
{
  if (domain) {
    //printf("domain %p & %p\n", *domain, domain);
    if (domain->destroy_folded)
      (*domain->destroy_folded)(domain->folded_params);
    if (domain->destroy_unfolded)
      (*domain->destroy_unfolded)(domain->unfolded_params);
    if (domain->destroy_k)
      (*domain->destroy_k)(domain->k_params);
    free(domain);
  }
}
@

<<destroy domain list>>=
void destroy_domain_list(list_t *domain_list)
{
  domain_list = head(domain_list);
  while (domain_list != NULL)
    destroy_domain((domain_t *) pop(&domain_list));
}
@

\subsection{Domain model and group handling}

<<group functions>>=
<<get tension handler>>
<<get group>>
<<int comparison function>>
<<find possible groups>>
<<is group active>>
<<get group list>>
<<get active groups>>
@ 

<<get tension handler>>=
one_dim_fn_t *get_tension_handler(domain_t *domain)
{
  if (domain->state == DS_FOLDED)
    return domain->folded_handler;
  else {
    if (domain->state != DS_UNFOLDED)
      fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
    assert(domain->state == DS_UNFOLDED);
    return domain->unfolded_handler;
  }
}
@ 

<<get group>>=
int get_group(domain_t *domain)
{
  if (domain->state == DS_FOLDED)
    return domain->folded_group;
  else {
    if (domain->state != DS_UNFOLDED)
      fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
    assert(domain->state == DS_UNFOLDED);
    return domain->unfolded_group;
  }
}
@ 

We already know all possible tension classes, since they are all
hardcoded into \prog.  However, there may be any number of possible
groups.  We define a function [[find_possible_groups]] to search for
possible (not neccessarily active) groups.  It can search for
subgroups of a single [[handler]], or by setting [[handler]] to
[[NULL]], subgroups of \emph{any} handler.  It returns a list of group
integers.
<<find possible groups>>=
list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
  list_t *ret = NULL;
  list = head(list);
  while (list != NULL) {
    if (handler == NULL || get_tension_handler(D(list)) == handler) {
      push(&ret, &D(list)->folded_group, 1);
      push(&ret, &D(list)->unfolded_group, 1);
    }
    list = list->next;
  }

  if (ret == NULL)
    return ret; /* no groups with this handler, no processing needed */

  /* sort the ret list, and remove duplicates */
  sort(&ret, &int_comparison);
  unique(&ret, &int_comparison);
#ifdef DEBUG
  fprintf(stderr, "handler %p possible groups:", handler);
  list = head(ret);
  while (list != NULL) {
    fprintf(stderr, " %d", *((int *)list->d));
    list = list->next;
  }
  fprintf(stderr, "\n");
#endif
  return ret;
}
@ 

Search a [[list]] of domains, and determine whether a particular model
class and group number combination is active.
<<is group active>>=
int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
{
  list = head(list);
  while (list != NULL) {
    if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
      return 1;
    list = list->next;
  }
  return 0;
}
@ 

Search a [[list]] of domains, and return all domains matching a
particular model class and group number.
<<get group list>>=
list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
{
  list_t *ret = NULL;
  list = head(list);
  while (list != NULL) {
    if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
      push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
#ifdef DEBUG
      fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
#endif
    }
    list = list->next;
  }
  return ret;
}
@ 
Because all the node data in lists returned by [[get_group_list]] is
also in the main domain list, you shouldn't destroy the node data
popped off when destroying the group lists.  It will all get cleaned
up when the main domain list is destroyed.

Put all this together to scan through a list of domains and construct
an array of all the active groups.
<<get active groups>>=
void get_active_groups(list_t *list,
                     int num_tension_handlers, one_dim_fn_t **tension_handlers,
                     int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
{
  void **active_groups=NULL;
  one_dim_fn_t *handler, **per_group_handlers=NULL;
  int i, num_active_groups, *group;
  list_t *possible_group_numbers,  *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;

  /* get all the active groups in a list */
  for (i=0; i<num_tension_handlers; i++) {
    handler = tension_handlers[i];
    if (handler == NULL) continue; /* tensionless handler */
    possible_group_numbers = head(find_possible_groups(list, handler));
    while (possible_group_numbers != NULL) {
      group = (int *)pop(&possible_group_numbers);
      if (is_group_active(list, handler, *group) == 1) {
        active_group_list = get_group_list(list, handler, *group);
        push(&active_groups_list, active_group_list, 1);
        push(&per_group_handlers_list, handler, 1);
#ifdef DEBUG
        fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
        list_print(stderr, active_group_list, "active group");
#endif
      }
    }
  }

  /* allocate the array we'll be returning */  
  num_active_groups = list_length(active_groups_list);
  active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
  assert(active_groups != NULL);
  per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
  assert(per_group_handlers != NULL);

  /* populate the array we'll be returning */
  active_groups_list = head(active_groups_list);
  for (i=0; i<num_active_groups; i++) {
    handler = pop(&per_group_handlers_list);
    per_group_handlers[i] = handler;

    active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
    assert(active_groups[i] != NULL);
    active_group_list = pop(&active_groups_list);
    ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
    ((tension_handler_data_t *)active_groups[i])->env = NULL;
    ((tension_handler_data_t *)active_groups[i])->persist = NULL;
  }
  assert(active_groups_list == NULL);
  assert(per_group_handlers_list == NULL);

  *pNum_active_groups = num_active_groups;
  *pPer_group_handlers = per_group_handlers;
  *pActive_groups = active_groups;
}
@ 


\section{String parsing}

For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
The model handling in getopt is set up to handle a fixed number of arguments for each model,
so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
need to take care of parsing those parameters themselves.
We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].

<<parse.h>>=
<<license comment>>
<<parse definitions>>
<<parse declarations>>
@ 

<<parse module makefile lines>>=
parse.c : sawsim.nw
        notangle -Rparse.c $^ > $@
parse.h : sawsim.nw
        notangle -Rparse.h $^ > $@
check_parse.c : sawsim.nw
        notangle -Rcheck-parse.c $^ > $@
check_parse : check_parse.c parse.c parse.h
        gcc -g -o $@ $< parse.c -lcheck
clean_parse : 
        rm -f parse.c parse.h check_parse.c check_parse
@ 

<<parse definitions>>=
#define SEP ',' /* argument separator character */
@ 

<<parse declarations>>=
extern void parse_list_string(char *string, char sep, char deeper, char shallower,
                       int *num, char ***string_array);
extern void free_parsed_list(int num, char **string_array);
@ 

[[parse_list_string]] allocates memory, don't forget to free it
afterward with [[free_parsed_list]].  It does not alter the original.

The string may start off with a [[deeper]] character
(i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
[[x,y]], where the pointer is one character in on the copied string.
However, when we go to free the memory, we need a pointer to the
beginning of the string.  In order to accommodate this for a string
with $N$ argument, allocate a pointer array with $N+1$ elements, let
the first $N$ elements point to the separated arguments, and let the
last element point to the start of the copied string regardless of
braces.
<<parse delimited list string functions>>=
/* TODO, split out into parse.hc */
static int next_delim_index(char *string, char sep, char deeper, char shallower)
{
  int i=0, depth = 0;
  while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
    if (string[i] == deeper) {depth++;}
    else if (string[i] == shallower) {depth--; assert(depth >= 0);}
    i++;
  }
  return i;
}

void parse_list_string(char *string, char sep, char deeper, char shallower,
                       int *num, char ***string_array)
{
  char *str=NULL, **ret=NULL;
  int i, j, n;
  if (string==NULL || strlen(string) == 0) {    /* handle the trivial cases */
    *num = 0;
    *string_array = NULL;
    return;
  }
  /* make a copy of the string, so we don't change the original */
  str = (char *)malloc(sizeof(char)*(strlen(string)+1));
  assert(str != NULL);
  strcpy(str, string); /* we know str is long enough */
  /* count the number of regions, so we can allocate pointers to them */
  i=-1; n=0;
  do {
    n++; i++; /* move on to next argument */
    i += next_delim_index(str+i, sep, deeper, shallower);
    //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
  } while (str[i] != '\0');
  ret = (char **)malloc(sizeof(char *)*(n+1));
  assert(ret != NULL);
  /* replace the separators with '\0' & assign pointers */
  ret[n] = str; /* point to the front of the copied string */
  j=0;
  ret[0] = str;
  for(i=1; i<n; i++) {
    j += next_delim_index(str+j, sep, deeper, shallower);
    str[j++] = '\0';
    ret[i] = str+j; /* point to the separated arguments */
  }
  /* strip off bounding braces, if any */
  for(i=0; i<n; i++) {
    if (ret[i][0] == deeper) {
      ret[i][0] = '\0';
      ret[i]++;
    }
    if (ret[i][strlen(ret[i])-1] == shallower)
      ret[i][strlen(ret[i])-1] = '\0';
  }
  *num = n;
  *string_array = ret;
}

void free_parsed_list(int num, char **string_array)
{
  if (string_array) {
    /* frees all items, since they were allocated as a single string*/
    free(string_array[num]);
    /* and free the array of pointers */
    free(string_array);
  }
}
@

\subsection{Parsing implementation}

<<parse.c>>=
<<license comment>>
<<parse includes>>
#include "parse.h"
<<parse delimited list string functions>>
@ 

<<parse includes>>=
#include <assert.h> /* assert()                */
#include <stdlib.h> /* NULL                    */
#include <stdio.h>  /* fprintf(), stdout       *//*!!*/
#include <string.h> /* strlen()                */
#include "parse.h"
@ 

\subsection{Parsing unit tests}

Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-parse.c>>=
<<parse check includes>>
<<string comparison definition and globals>>
<<check relative error>>
<<parse test suite>>
<<main check program>>
@ 

<<parse check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
#include <stdio.h>  /* printf()                              */
#include <assert.h> /* assert()                              */
#include <string.h> /* strlen()                              */
<<check includes>>
#include "parse.h"
@ 

<<parse test suite>>=
<<parse list string tests>>
<<parse suite definition>>
@ 

<<parse suite definition>>=
Suite *test_suite (void)
{
  Suite *s = suite_create ("k model");
  <<parse list string test case defs>>

  <<parse list string test case add>>
  return s;
}
@ 

<<parse list string tests>>=

/*
START_TEST(test_next_delim_index)
{
  fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
  fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
  fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
  fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
  fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
}
END_TEST
*/
START_TEST(test_parse_list_null)
{
  int num_param_args;
  char **param_args;
  parse_list_string(NULL, SEP, '{', '}',
                    &num_param_args, &param_args);
  fail_unless(num_param_args == 0, NULL);
  fail_unless(param_args == NULL, NULL);
}
END_TEST
START_TEST(test_parse_list_single_simple)
{
  int num_param_args;
  char **param_args;
  parse_list_string("arg", SEP, '{', '}',
                    &num_param_args, &param_args);
  fail_unless(num_param_args == 1, NULL);
  fail_unless(STRMATCH(param_args[0],"arg"), NULL);
}
END_TEST
START_TEST(test_parse_list_single_compound)
{
  int num_param_args;
  char **param_args;
  parse_list_string("{x,y,z}", SEP, '{', '}',
                    &num_param_args, &param_args);
  fail_unless(num_param_args == 1, NULL);
  fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
}
END_TEST
START_TEST(test_parse_list_double_simple)
{
  int num_param_args;
  char **param_args;
  parse_list_string("abc,def", SEP, '{', '}',
                    &num_param_args, &param_args);
  fail_unless(num_param_args == 2, NULL);
  fail_unless(STRMATCH(param_args[0],"abc"), NULL);
  fail_unless(STRMATCH(param_args[1],"def"), NULL);
}
END_TEST
@
<<parse list string test case defs>>=
TCase *tc_parse_list_string = tcase_create("parse list string");
@
<<parse list string test case add>>=
//tcase_add_test(tc_parse_list_string, test_next_delim_index);
tcase_add_test(tc_parse_list_string, test_parse_list_null);
tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
suite_add_tcase(s, tc_parse_list_string);
@


\section{Unit tests}

Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<checks>>=
<<includes>>
<<check includes>>
<<definitions>>
<<globals>>
<<check globals>>
<<functions>>
<<test suite>>
<<main check program>>
@ 

<<check includes>>=
#include <check.h>
@ 

<<check globals>>=
@ 

<<test suite>>=
<<F tests>>
<<determine dt tests>>
<<happens tests>>
<<does domain unfold tests>>
<<randomly unfold domains tests>>
<<suite definition>>
@ 

<<suite definition>>=
Suite *test_suite (void)
{
  Suite *s = suite_create ("sawsim");
  <<F test case defs>>
  <<determine dt test case defs>>
  <<happens test case defs>>
  <<does domain unfold test case defs>>
  <<randomly unfold domains test case defs>>

  <<F test case add>>
  <<determine dt test case add>>
  <<happens test case add>>
  <<does domain unfold test case add>>
  <<randomly unfold domains test case add>>

/*
  tcase_add_checked_fixture(tc_strip_address,
                            setup_strip_address,
                            teardown_strip_address);
*/
  return s;
}
@ 

<<main check program>>=
int main(void)
{
  int number_failed;
  Suite *s = test_suite();
  SRunner *sr = srunner_create(s);
  srunner_run_all(sr, CK_ENV);
  number_failed = srunner_ntests_failed(sr);
  srunner_free(sr);
  return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}
@ 

\subsection{F tests}
<<F tests>>=
<<worm-like chain tests>>
@ 
<<F test case defs>>=
<<worm-like chain test case def>>
@ 
<<F test case add>>=
<<worm-like chain test case add>>
@ 

<<worm-like chain tests>>=
START_TEST(test_wlc_at_zero)
{
  double T=1.0, L=1.0, p=0.1, x=0.0;
  fail_unless(wlc(x, T, p, L)==0, NULL);
}
END_TEST

START_TEST(test_wlc_at_half)
{
  double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
  /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
   * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
   *                = 0.25 * 3 = 3/4
   * linear term = x/L = 0.5
   * nonlinear + linear = 0.75 + 0.5 = 1.25
   * wlc = 10e21*1.25 = 12.5e21 
   */
  fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
              "wlc(%g, %g, %g, %g) = %g != %g",
               x, T, p, L, wlc(x, T, p, L), 12.5e21);
}
END_TEST
@ 
<<worm-like chain test case def>>=
TCase *tc_wlc = tcase_create("WLC");
@ 

<<worm-like chain test case add>>=
tcase_add_test(tc_wlc, test_wlc_at_zero);
tcase_add_test(tc_wlc, test_wlc_at_half);
suite_add_tcase(s, tc_wlc);
@ 

\subsection{Model tests}

Check the searching with [[linear_k]].
Check overwhelming force treatment with the heavyside-esque step [[?]].
<<determine dt tests>>=
double linear_k(double F, environment_t *env, void *params)
{
  double Fnot = *(double *)params;
  return Fnot+F;
}

START_TEST(test_determine_dt_linear_k)
{
  environment_t env;
  double dt_max=3.0, Fnot=3.0;
  double F[]={0,1,2,3,4,5,6};
  domain_t dom; /* use both parts at once for folded/unfolded */
  int i;
  env.T = 300.0;
/*
  dom->next = dom->prev = NULL;
  dom->k_func_t = linear_k;
  dom->folded_params = &Fnot;
  dom->unfolded_params = !!!!!!!!!
  dom->destroy_folded = dom->destroy_unfolded = NULL;
  for( i=0; i < sizeof(F)/sizeof(double); i++) {
    fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
  }
*/
}
END_TEST
@
<<determine dt test case defs>>=
TCase *tc_determine_dt = tcase_create("Determine dt");
@
<<determine dt test case add>>=
tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
suite_add_tcase(s, tc_determine_dt);
@

<<happens tests>>=
@
<<happens test case defs>>=
@
<<happens test case add>>=
@

<<does domain unfold tests>>=
@
<<does domain unfold test case defs>>=
@
<<does domain unfold test case add>>=
@

<<randomly unfold domains tests>>=
@
<<randomly unfold domains test case defs>>=
@
<<randomly unfold domains test case add>>=
@


\section{Balancing group extensions}
\label{app.tension_balance}

For a given total extension $x$ (of the piezo), the various domain
models (WLC, FJC, Hookean springs, \ldots) and groups extend different
amounts, and we need to tweak the portion that each extends to balance
the tension amongst the active groups.  Since the tension balancing is
separable from the bulk of the simulation, we break it out into a
separate module.  The interface is defined in [[tension_balance.h]],
the implementation in [[tension_balance.c]], and the unit testing in
[[check_tension_balance.c]]

<<tension-balance.h>>=
<<license comment>>
<<tension balance function declaration>>
@ 

<<tension balance functions>>=
<<one dimensional bracket>>
<<one dimensional solve>>
<<x of x0>>
<<tension balance function>>
@ 

<<tension balance module makefile lines>>=
tension_balance.c : sawsim.nw
        notangle -Rtension-balance.c $^ > $@
tension_balance.h : sawsim.nw
        notangle -Rtension-balance.h $^ > $@
check_tension_balance.c : sawsim.nw
        notangle -Rcheck-tension-balance.c $^ > $@
check_tension_balance : check_tension_balance.c global.h tension_balance.c tension_balance.h
        gcc -g -o $@ $< tension_balance.c -lcheck
clean_tension : 
        rm -f tension_balance.c tension_balance.h check_tension_balance.c check_tension_balance
@ 

The entire force balancing problem reduces to a solving two nested
one-dimensional root-finding problems.  First we define the one
dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
non-empty group).  $x(x_0)$ is also strictly monotonically increasing,
so we can solve for $x_0$ such that $\sum_i x_i = x$.  [[last_x]]
stores the last successful value of $x$, since we expect to be taking
small steps ($x \approx$ [[last_x]]).  Setting [[last_x = -1]]
indicates that the groups have changed.
<<tension balance function declaration>>=
double tension_balance(int num_tension_groups,
                       one_dim_fn_t **tension_handlers,
                       void **params,
                       double last_x,
                       double x);
@ 
<<tension balance function>>=
double tension_balance(int num_tension_groups,
                       one_dim_fn_t **tension_handlers,
                       void **params,
                       double last_x,
                       double x)
{
  static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
  double F, xo_guess, xo, lb, ub;
  double min_dx=1e-15, min_dy=1e-16;
  int max_steps=100, i;

  assert(num_tension_groups > 0);
  assert(tension_handlers != NULL);
  assert(params != NULL);
  assert(num_tension_groups > 0);

  if (num_tension_groups == 1) { /* only one group, no balancing required */
    xo = x;
  } else {
    //fprintf(stderr, "balancing for x = %g with ", x);
    //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
    //fprintf(stderr, "\n");
    if (last_x == -1) { /* new group setup, reset x_xo_data */
      /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
      if (x_xo_data.xi != NULL)
        free(x_xo_data.xi);
      /* construct new x_xo_data */
      x_xo_data.n_groups = num_tension_groups;
      x_xo_data.tension_handler = tension_handlers;
      x_xo_data.handler_data = params;
      x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
      for (i=0; i<num_tension_groups; i++)
        x_xo_data.xi[i] = -1.0;
    }
    if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
      if (x == 0) xo_guess = length_scale;
      else        xo_guess = x/num_tension_groups;
#ifdef DEBUG
      fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
#endif
      oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
    } else { /* work off the last balanced point */
#ifdef DEBUG
      fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
#endif
      if (x > last_x) {
        lb = x_xo_data.xi[0];
        ub = x_xo_data.xi[0]+ x-last_x;   /* apply all change to x0 */
      } else if (x < last_x) {
        lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
        ub = x_xo_data.xi[0];
      } else { /* x == last_x */
        fprintf(stderr, "not moving\n");
        lb= x_xo_data.xi[0];
        ub= x_xo_data.xi[0];
      }
    }
#ifdef DEBUG
    fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
            __FUNCTION__, x, lb, ub);
#endif
    xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_dx, min_dy, max_steps, NULL);
#ifdef DEBUG
    if (num_tension_groups > 1) {
      fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
      for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
      fprintf(stderr, "\n");
    }
#endif
  }
  F = (*tension_handlers[0])(xo, params[0]);
  return F;
}
@ 

<<tension balance internal definitions>>=
<<x of x0 definitions>>
@ 

<<x of x0 definitions>>=
typedef struct x_of_xo_data_struct {
  int n_groups;
  one_dim_fn_t **tension_handler; /* array of fn pointers */
  void **handler_data;            /* array of void* pointers */
  double *xi;
} x_of_xo_data_t;
@ 

<<x of x0>>=
double x_of_xo(double xo, void *pdata)
{
  x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
  double F, x=0, xi, xi_guess, lb, ub;
  int i;
  double min_dx=1e-15, min_dy=1e-16;
  int max_steps=100;
  assert(data->n_groups > 0);
  data->xi[0] = xo;
  F = (*data->tension_handler[0])(xo, data->handler_data[0]);
  x += xo;
  if (data->xi)  data->xi[0] = xo;
  for (i=1; i < data->n_groups; i++) {
    if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
    else                   xi_guess = data->xi[i];
    oneD_bracket(data->tension_handler[i], data->handler_data[i], F,  xi_guess, &lb, &ub);
    xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_dx, min_dy, max_steps, NULL);
    data->xi[i] = xi;
    x += xi;
    if (data->xi)  data->xi[i] = xi;
  }
  return x;
}
@ 

Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited number of steps for monotonically increasing $f(x)$.
Simple bisection, so it's robust and converges fairly quickly.
<<one dimensional function definition>>=
/* equivalent to gsl_func_t */
typedef double one_dim_fn_t(double x, void *params);
@ 
<<one dimensional solve declaration>>=
double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
                  double min_dx, double min_dy, int max_steps, 
                  double *pYx);
@ 

<<one dimensional solve>>=
/* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
                  double min_dx, double min_dy, int max_steps, 
                  double *pYx)
{
  double yx, ylb, yub, x;
  int i=0;
  assert(ub >= lb);
  ylb = (*f)(lb, params);
  yub = (*f)(ub, params);
  
  /* check some simple cases */
  if (ylb == yub) {
    if (ylb != y)  return HUGE_VAL; /* error! f(x)=y not bounded */
    else           return lb; /* any x on the range [lb, ub] would work */
  }
  if (ylb == y) { x=lb; yx=ylb; goto end; }
  if (yub == y) { x=ub; yx=yub; goto end; }

#ifdef DEBUG
  fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
#endif
  assert(ylb < y); assert(yub > y);
  for (i=0; i<max_steps; i++) {
    x = (lb+ub)/2.0;
    yx = (*f)(x, params);
#ifdef DEBUG
  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);
#endif
    if (ub-lb < min_dx || yub-ylb < min_dy || yx == y) break;
    else if (yx > y)  { ub=x; yub=yx; }
    else /* yx < y */ { lb=x; ylb=yx; }
  }
 end:
  if (pYx) *pYx = yx;
#ifdef DEBUG
  fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
#endif
  return x;
}
@ 

The one dimensional solver needs a bracketed solution, which sometimes we don't have.
Generate bracketing $x$ values through bisection or exponential growth.
<<one dimensional bracket declaration>>=
void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
@ 

<<one dimensional bracket>>=
void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
{
  double yx, step, x=xguess;
  int i=0;
  yx = (*f)(x, params);
#ifdef DEBUG
  fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
#endif
  if (yx > y) {
    assert(x > 0.0);
    *ub = x;
    *lb = 0;
#ifdef DEBUG
    fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
#endif
  } else {
    *lb = x;
    if (x == 0) x = length_scale/2.0; /* HACK */
    while (yx < y) {
      x *= 2.0;
      yx = (*f)(x, params);
#ifdef DEBUG
      fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
#endif
    }
    *ub = x;
#ifdef DEBUG
    fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
#endif
  }
  //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
}
@ 

\subsection{Balancing implementation}

<<tension-balance.c>>=
//#define DEBUG
<<license comment>>
<<tension balance includes>>
#include "tension_balance.h"

double length_scale = 1e-10; /* HACK */

<<tension balance internal definitions>>
<<tension balance functions>>
@ 

<<tension balance includes>>=
#include <assert.h> /* assert()          */
#include <stdlib.h> /* NULL              */
#include <math.h>   /* HUGE_VAL, macro constant, so don't need to link to math */
#include <stdio.h>  /* fprintf(), stdout */
#include "global.h"
@ 

\subsection{Balancing unit tests}

Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-tension-balance.c>>=
<<tension balance check includes>>
<<tension balance test suite>>
<<main check program>>
@ 

<<tension balance check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
#include <assert.h> /* assert()                      */
<<check includes>>
#include "global.h"
#include "tension_balance.h"
@ 

<<tension balance test suite>>=
<<tension balance function tests>>
<<tension balance suite definition>>
@ 

<<tension balance suite definition>>=
Suite *test_suite (void)
{
  Suite *s = suite_create ("tension balance");
  <<tension balance function test case defs>>

  <<tension balance function test case adds>>
  return s;
}
@ 

<<tension balance function tests>>=
<<check relative error>>

double hooke(void *pK, double x)
{
  assert(x >= 0);
  return *((double*)pK) * x;
}

START_TEST(test_single_function)
{
  double k=5, x=3, last_x=2.0, F;
  one_dim_fn_t *handlers[] = {&hooke};
  void *data[] =               {&k};
  double xi[] =                {0.0};
  int active[] =               {1};
  int new_active_groups = 1;
  F = tension_balance(1, handlers, data, active, new_active_groups, xi, last_x, x);
  fail_unless(F = k*x, NULL);
}
END_TEST
@ 

We can also test balancing two springs with different spring constants.
The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
<<tension balance function tests>>=
START_TEST(test_double_hooke)
{
  double k1=5, k2=4, x=3, last_x=2.0, F, Fe, x1e, x2e;
  one_dim_fn_t *handlers[] = {&hooke, &hooke, NULL};
  void *data[] =               {&k1,    &k2,    NULL};
  double xi[] =                {0.0,    0.0,    0.0};
  int active[] =               {1,      1,      0};
  int new_active_groups = 1;
  F = tension_balance(3, handlers, data, active, new_active_groups, xi, last_x, x);
  x1e = x*k2/(k1+k2);
  Fe = k1*x1e;
  x2e = x1e*k1/k2;
  //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
  CHECK_ERR(1e-6, x1e, xi[0]);
  CHECK_ERR(1e-6, x2e, xi[1]);
  CHECK_ERR(1e-6, Fe, F);
}
END_TEST
@ 

<<tension balance function test case defs>>=
TCase *tc_tbfunc = tcase_create("tension balance function");
@ 

<<tension balance function test case adds>>=
tcase_add_test(tc_tbfunc, test_single_function);
tcase_add_test(tc_tbfunc, test_double_hooke);
suite_add_tcase(s, tc_tbfunc);
@ 

\section{Lists}

The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]

<<list.h>>=
<<license comment>>
<<list definitions>>
<<list declarations>>
@ 

<<list declarations>>=
<<head and tail declarations>>
<<list length declaration>>
<<push declaration>>
<<pop declaration>>
<<list destroy declaration>>
<<list to array declaration>>
<<move declaration>>
<<sort declaration>>
<<unique declaration>>
<<list print declaration>>
@ 

<<list functions>>=
<<create and destroy node>>
<<head and tail>>
<<list length>>
<<push>>
<<pop>>
<<list destroy>>
<<list to array>>
<<move>>
<<sort>>
<<unique>>
<<list print>>
@ 

<<list module makefile lines>>=
list.c : sawsim.nw
        notangle -Rlist.c $^ > $@
list.h : sawsim.nw
        notangle -Rlist.h $^ > $@
check_list.c : sawsim.nw
        notangle -Rcheck-list.c $^ > $@
check_list : check_list.c global.h list.c list.h
        gcc -g -o $@ $< list.c -lcheck
clean_list : 
        rm -f list.c list.h check_list.c check_list
@ 

<<list definitions>>=
typedef struct list_struct {
  struct list_struct *next;
  struct list_struct *prev;
  void *d; /* data */
} list_t;

<<comparison function typedef>>
@ 

[[head]] and [[tail]] return pointers to the head and tail nodes of the list:
<<head and tail declarations>>=
list_t *head(list_t *list);
list_t *tail(list_t *list);
@ 
<<head and tail>>=
list_t *head(list_t *list)
{
  if (list == NULL)
    return list;
  while (list->prev != NULL) {
    list = list->prev;
  }
  return list;
}

list_t *tail(list_t *list)
{
  if (list == NULL)
    return list;
  while (list->next != NULL) {
    list = list->next;
  }
  return list;
}
@ 

<<list length declaration>>=
int list_length(list_t *list);
@ 
<<list length>>=
int list_length(list_t *list)
{
  int i;
  if (list == NULL)
    return 0;
  list = head(list);
  i = 1;
  while (list->next != NULL) {
    list = list->next;
    i += 1;
  }
  return i;
}
@ 

[[push]] inserts a node relative to the current position in the list
without changing the current position.
However, if the list we're pushing onto is [[NULL]], the current position isn't defined
so we set it to that of the pushed domain.
If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
<<push declaration>>=
void push(list_t **pList, void *data, int below);
@ 
<<push>>=
void push(list_t **pList, void *data, int below)
{  
  list_t *list, *new_node;
  assert(pList != NULL);
  list = *pList;
  new_node = create_node(data);
  if (list == NULL)
    *pList = new_node;
  else if (below==0) { /* insert above */
    if (list->prev != NULL)
      list->prev->next = new_node;
    new_node->prev = list->prev;
    list->prev = new_node;
    new_node->next = list;
  } else {         /* insert below */
    if (list->next != NULL)
      list->next->prev = new_node;
    new_node->next = list->next;
    list->next = new_node;
    new_node->prev = list;
  }
}
@ 

[[pop]] removes the current domain node, moving the current position
to the node after it, unless that node is [[NULL]], in which case move
the current position to the node before it.
<<pop declaration>>=
void *pop(list_t **pList);
@ 
<<pop>>=
void *pop(list_t **pList)
{
  list_t *list, *popped;
  void *data;
  assert(pList != NULL);
  list = *pList;
  assert(list != NULL); /* not an empty list */
  popped = list;
  /* bypass the popped node */
  if (list->prev != NULL)
     list->prev->next = list->next;
  if (list->next != NULL) {
     list->next->prev = list->prev;
     *pList = list->next;
  } else
     *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
  data = popped->d;
  destroy_node(popped);
  return data;
}
@ 

[[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
If the cleanup function is [[NULL]], no cleanup function is called.
The nodes are not popped in any particular order.
<<list destroy declaration>>=
void list_destroy(list_t **pList, destroy_data_func_t *destroy);
@ 
<<list destroy>>=
void list_destroy(list_t **pList, destroy_data_func_t *destroy)
{
  list_t *list;
  void *data;
  assert(pList != NULL);
  list = *pList;
  *pList = NULL;
  assert(list != NULL); /* not an empty list */
  while (list != NULL) {
    data = pop(&list);
    if (destroy != NULL)
      destroy(data);
  }
}
@ 

[[list_to_array]] converts a list to an array of pointers, preserving order.
<<list to array declaration>>=
void list_to_array(list_t *list, int *length, void ***array);
@ 
<<list to array>>=
void list_to_array(list_t *list, int *pLength, void ***pArray)
{
  void **array;
  int i,length;
  assert(list != NULL);
  assert(pLength != NULL);
  assert(pArray != NULL);

  length = list_length(list);
  array = (void **)malloc(sizeof(void *)*length);
  assert(array != NULL);
  list = head(list);
  i=0;
  while (list != NULL) {
    array[i++] = list->d;
    list = list->next;
  }
  *pLength = length;
  *pArray = array;
}
@ 

[[move]] moves the current element along the list in either direction.
<<move declaration>>=
void move(list_t **pList, int downstream);
@ 
<<move>>=
void move(list_t **pList, int downstream)
{
  list_t *list, *flipper;
  void *data;
  assert(pList != NULL);
  list = *pList;          /* a>B>c>d */
  assert(list != NULL); /* not an empty list */
  if (downstream == 0)
    flipper = list->prev; /* flipper is A */
  else
    flipper = list->next; /* flipper is C */
  /* ensure there is somewhere to go in the direction we're heading */
  assert(flipper != NULL);
  /* remove the current element */
  data = pop(&list);      /* data=B, a>C>d */
  /* and put it back in in the correct spot */
  list = flipper;
  if (downstream == 0) {
    push(&list, data, 0); /* b>A>c>d */
    list = list->prev;    /* B>a>c>d   */    
  } else {
    push(&list, data, 1); /* a>C>b>d */
    list = list->next;    /* a>c>B>d */
  }
  *pList = list;
}
@ 

[[sort]] sorts a list from smallest to largest according to a user
specified comparison.
<<comparison function typedef>>=
typedef int (comparison_fn_t)(void *A, void *B);
@ 
For example
<<int comparison function>>=
int int_comparison(void *A, void *B)
{
  int a,b;
  a = *((int *)A);
  b = *((int *)B);
  if (a > b) return 1;
  else if (a == b) return 0;
  else return -1;
}
@ 
<<sort declaration>>=
void sort(list_t **pList, comparison_fn_t *comp);
@ 
Here we use the bubble sort algorithm.
<<sort>>=
void sort(list_t **pList, comparison_fn_t *comp)
{
  list_t *list;
  int stable=0;
  assert(pList != NULL);
  list = *pList;
  assert(list != NULL); /* not an empty list */
  while (stable == 0) {
    stable = 1;
    list = head(list);
    while (list->next != NULL) {
      if ((*comp)(list->d, list->next->d) > 0) {
        stable = 0;
        move(&list, 1 /* downstream */);
      } else
        list = list->next;
    }
  }
  *pList = list;
}
@ 

[[unique]] removes duplicates from a sorted list according to a user
specified comparison.  Don't do this unless you have other pointers
any dynamically allocated data in your list, because [[unique]] just
drops any non-unique data without freeing it.
<<unique declaration>>=
void unique(list_t **pList, comparison_fn_t *comp);
@ 
<<unique>>=
void unique(list_t **pList, comparison_fn_t *comp)
{
  list_t *list;
  assert(pList != NULL);
  list = *pList;
  assert(list != NULL); /* not an empty list */
  list = head(list);
  while (list->next != NULL) {
    if ((*comp)(list->d, list->next->d) == 0) {
      pop(&list);
    } else
      list = list->next;
  }
  *pList = list;
}
@ 

[[list_print]] prints a list to a [[FILE *]] stream.
<<list print declaration>>=
void list_print(FILE *file, list_t *list, char *list_name);
@ 
<<list print>>=
void list_print(FILE *file, list_t *list, char *list_name)
{
  fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
  list = head(list);
  while (list != NULL) {
    fprintf(file, " %p", list->d);
    list = list->next;
  }
  fprintf(file, "\n");
}
@ 

[[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
<<create and destroy node>>=
list_t *create_node(void *data)
{
  list_t *ret = (list_t *)malloc(sizeof(list_t));
  assert(ret != NULL);
  ret->prev = NULL;
  ret->next = NULL;
  ret->d = data;
  return ret;
}

void destroy_node(list_t *node)
{
  if (node)
    free(node);
}
@
The user must free the data pointed to by the node on their own.

\subsection{List implementation}

<<list.c>>=
<<license comment>>
<<list includes>>
#include "list.h"
<<list functions>>
@ 

<<list includes>>=
#include <assert.h> /* assert()                                */
#include <stdlib.h> /* malloc(), free(), rand()                */
#include <stdio.h>  /* fprintf(), stdout, FILE                 */
#include "global.h" /* destroy_data_func_t */
@ 

\subsection{List unit tests}

Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-list.c>>=
<<list check includes>>
<<list test suite>>
<<main check program>>
@ 

<<list check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
<<check includes>>
#include "list.h"
@ 

<<list test suite>>=
<<push tests>>
<<pop tests>>
<<list suite definition>>
@ 

<<list suite definition>>=
Suite *test_suite (void)
{
  Suite *s = suite_create ("list");
  <<push test case defs>>
  <<pop test case defs>>

  <<push test case adds>>
  <<pop test case adds>>
  return s;
}
@ 

<<push tests>>=
START_TEST(test_push)
{
  list_t *list=NULL;
  int i, p, e, n=3;
  for (i=0; i<n; i++)
    push(&list, (void *)i, 1);
  fail_unless(list == head(list), NULL);
  fail_unless( (int)list->d == 0 );
  for (i=0; i<n; i++) {
    /* we expect 0, n-1, n-2, ..., 1 */
    if (i == 0) e = 0;
    else        e = n-i;
    fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
  }
}
END_TEST
@ 

<<push test case defs>>=
TCase *tc_push = tcase_create("push");
@ 

<<push test case adds>>=
tcase_add_test(tc_push, test_push);
suite_add_tcase(s, tc_push);
@ 

<<pop tests>>=
@ 
<<pop test case defs>>=
@ 
<<pop test case adds>>=
@ 

\section{Function string evaluation}

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).
We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
The solution is to run a scripting language as a subprocess accessed via pipes.
We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.

We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
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.
We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.

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.
Then you can either statically or dynamically link to those hardcoded functions.
While much less flexible, this approach would be a fairly simple-to-implement fallback.

Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].

<<string-eval.h>>=
<<license comment>>
<<string eval setup declaration>>
<<string eval function declaration>>
<<string eval teardown declaration>>
@ 

<<string eval module makefile lines>>=
string_eval.c : sawsim.nw
        notangle -Rstring-eval.c $^ > $@
string_eval.h : sawsim.nw
        notangle -Rstring-eval.h $^ > $@
check_string_eval.c : sawsim.nw
        notangle -Rcheck-string-eval.c $^ > $@
check_string_eval : check_string_eval.c string_eval.c string_eval.h
        gcc -g -o $@ $< string_eval.c -lcheck -lgsl -lgslcblas -lm
clean_string_eval :
        rm -f string_eval.c string_eval.h check_string_eval.c check_string_eval
@ 

For an introduction to POSIX process control, see\\
 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
 and of course, the relavent [[man]] pages.

We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
[[execvp]] replaces the calling process' program with a new program.
The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
(this may be a security hole if someone messes about with [[PATH]] before you run \prog,
 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
The new program needs command line arguments, just like it would if you were running it from a shell.
The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
with the final array entry being a [[NULL]] pointer.

Now that we know how [[execvp]] works, we store it's arguments in some definitions,
to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
<<string eval subprocess definitions>>=
#define SUBPROCESS "python"
//#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
//static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
@ 
The [[i]] option lets Python know that it should run in interactive mode.
In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
In interactive mode, python acts on every instruction as soon as it is recieved.
The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply turns off the default prompting ([[ps1]] and [[ps2]]), since that is mostly for human users, and our program doesn't need it. 
%The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply [[pass]]es so that the silly Python header information is not printed.  We leave the prompts in, because we scan for them to determine when the output has completed.

Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
[[fork]] returns two copies of the same program, executing the original code.
In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].

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.
The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
<<string eval pipe definitions>>=
#define PIPE_READ  0   /* the end you read from */
#define PIPE_WRITE 1   /* the end you write to */

#define STDIN      0   /* initial index of stdin pair  */
#define STDOUT     2   /* initial index of stdout pair */

@ 
So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.

As a finishing touch, we can promote the POSIX file descriptors ([[read]]/[[write]] interface) into the more familiar [[stdio.h]] \emph{streams} ([[fprintf]]/[[fgetc]] interface) using [[fdopen]], which creates a stream from an open file descriptor.

<<string eval setup declaration>>=
extern void string_eval_setup(FILE **pIn, FILE **pOut);
@ 
<<string eval setup definition>>=
void string_eval_setup(FILE **pIn, FILE **pOut)
{
  pid_t pid;
  int pfd[4]; /* file descriptors for stdin and stdout */
  int rc;
  assert(pipe(pfd+STDIN)  != -1); /* stdin pair  (in, out) */
  assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
  
  pid = fork(); /* split process into two copies */

  if (pid == -1) { /* fork error */
    perror("fork error");
    exit(1);
  } else if (pid == 0) { /* child */
    close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input   */
    close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
    dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin  (closes original stdin) */
    dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
    execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
    perror("exec error"); /* exec shouldn't return */
    _exit(1);
  } else { /* parent */
    close(pfd[STDIN+PIPE_READ]);   /* close stdin pipe output */
    close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
    *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
    if ( *pIn == NULL ) {
      perror("fdopen (in)");
      exit(1);
    }
    *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
    if ( *pOut == NULL ) {
      perror("fdopen (out)");
      exit(1);
    }
  }
}
@

To use the evaluating subprocess, we just pipe in our command, and read out the results.
For the simple cases we expect here, we restrict ourselves to a single line of returned text.
<<string eval function declaration>>=
extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
@
<<string eval function definition>>=
void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
{
  int rc;
  rc = fprintf(in, "%s", input);
  assert(rc == strlen(input));
  fflush(in);
  fflush(out);
  alarm(1); /* set a one second timeout on the read */
  assert( fgets(output, buflen, out) != NULL );
  alarm(0); /* cancel the timeout */
  //fprintf(stderr, "eval: %s ----> %s", input, output);
}
@
The [[alarm]] calls set and clear a timeout on the returned output.
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.
This protects against invalid input for which a line of output is not printed to [[stdout]].
Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
If you are getting strange results, check your python code seperately. TODO, better error handling.

Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
<<string eval teardown declaration>>=
extern void string_eval_teardown(FILE **pIn, FILE **pOut);
@ 

<<string eval teardown definition>>=
void string_eval_teardown(FILE **pIn, FILE **pOut)
{
  pid_t pid=0;
  int stat_loc;

  /* redirect Python's stderr */
  fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
  fflush(*pIn);

  /* close pipes */
  assert( fclose(*pIn) == 0 );
  *pIn = NULL;
  assert( fclose(*pOut) == 0 );
  *pOut = NULL;

  /* wait for python to exit */
  while (pid <= 0) {
    pid = wait(&stat_loc);
    if (pid < 0) {
      perror("pid");
    }
  }

  /*
  if (WIFEXITED(stat_loc)) {
    printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
  } else if (WIFSIGNALED(stat_loc)) {
    printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
  }
  */
}
@ 
The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.

\subsection{String evaluation implementation}

<<string-eval.c>>=
<<license comment>>
<<string eval includes>>
#include "string_eval.h"
<<string eval internal definitions>>
<<string eval functions>>
@ 

<<string eval includes>>=
#include <assert.h> /* assert()                    */
#include <stdlib.h> /* NULL                        */
#include <stdio.h>  /* fprintf(), stdout, fdopen() */
#include <string.h> /* strlen()                    */
#include <sys/types.h> /* pid_t                    */
#include <unistd.h> /* pipe(), fork(), execvp(), alarm()    */
#include <wait.h>   /* wait()                      */
@ 

<<string eval internal definitions>>=
<<string eval subprocess definitions>>
<<string eval pipe definitions>>
@ 

<<string eval functions>>=
<<string eval setup definition>>
<<string eval function definition>>
<<string eval teardown definition>>
@ 

\subsection{String evaluation unit tests}

Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-string-eval.c>>=
<<string eval check includes>>
<<string comparison definition and globals>>
<<string eval test suite>>
<<main check program>>
@ 

<<string eval check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
#include <stdio.h>  /* printf()                              */
#include <assert.h> /* assert()                              */
#include <string.h> /* strlen()                              */
#include <signal.h> /* SIGKILL                               */
<<check includes>>
#include "string_eval.h"
@ 

<<string eval test suite>>=
<<string eval tests>>
<<string eval suite definition>>
@ 

<<string eval suite definition>>=
Suite *test_suite (void)
{
  Suite *s = suite_create ("string eval");
  <<string eval test case defs>>

  <<string eval test case add>>
  return s;
}
@ 

<<string eval tests>>=
START_TEST(test_setup_teardown)
{
  FILE *in, *out;
  string_eval_setup(&in, &out);
  string_eval_teardown(&in, &out);
}
END_TEST
START_TEST(test_invalid_command)
{
  FILE *in, *out;
  char input[80], output[80]={};
  string_eval_setup(&in, &out);
  sprintf(input, "print ABCDefg\n");
  string_eval(in, out, input, 80, output); 
  string_eval_teardown(&in, &out);
}
END_TEST
START_TEST(test_simple_eval)
{
  FILE *in, *out;
  char input[80], output[80]={};
  string_eval_setup(&in, &out);
  sprintf(input, "print 3+4*5\n");
  string_eval(in, out, input, 80, output); 
  fail_unless(STRMATCH(output,"23\n"), NULL);
  string_eval_teardown(&in, &out);
}
END_TEST
START_TEST(test_multiple_evals)
{
  FILE *in, *out;
  char input[80], output[80]={};
  string_eval_setup(&in, &out);
  sprintf(input, "print 3+4*5\n");
  string_eval(in, out, input, 80, output); 
  fail_unless(STRMATCH(output,"23\n"), NULL);
  sprintf(input, "print (3**2 + 4**2)**0.5\n");
  string_eval(in, out, input, 80, output); 
  fail_unless(STRMATCH(output,"5.0\n"), NULL);
  string_eval_teardown(&in, &out);
}
END_TEST
START_TEST(test_eval_with_state)
{
  FILE *in, *out;
  char input[80], output[80]={};
  string_eval_setup(&in, &out);
  sprintf(input, "print 3+4*5\n");
  fprintf(in, "A = 3\n");
  sprintf(input, "print A*3\n");
  string_eval(in, out, input, 80, output); 
  fail_unless(STRMATCH(output,"9\n"), NULL);
  string_eval_teardown(&in, &out);
}
END_TEST
@
<<string eval test case defs>>=
TCase *tc_string_eval = tcase_create("string_eval");
@
<<string eval test case add>>=
tcase_add_test(tc_string_eval, test_setup_teardown);
tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
tcase_add_test(tc_string_eval, test_simple_eval);
tcase_add_test(tc_string_eval, test_multiple_evals);
tcase_add_test(tc_string_eval, test_eval_with_state);
suite_add_tcase(s, tc_string_eval);
@


\section{Accelerating function evaluation}

My first version-0.3 code was running very slowly.
With the options suggested in the help 
([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]), 
the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
That's only 15 calls per solution, so the search algorithm seems reasonable.
The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.

<<accel-k.h>>=
double accel_k(k_func_t *k, double F, environment_t *env, void *params);
void free_accels();
@ 

<<accel k module makefile lines>>=
accel_k.c : sawsim.nw
        notangle -Raccel-k.c $^ > $@
accel_k.h : sawsim.nw
        notangle -Raccel-k.h $^ > $@
check_accel_k.c : sawsim.nw
        notangle -Rcheck-accel_k.c $^ > $@
check_accel_k : check_accel_k.c global.h
        gcc -g -o $@ $< accel_k.c -lcheck -lgsl -lgslcblas -lm
clean_accel_k : 
        rm -f accel_k.c accel_k.h check_accel_k.c check_accel_k
@ 

<<accel-k.c>>=
#include <assert.h>  /* assert()                */
#include <stdlib.h>  /* realloc(), free(), NULL */
#include "global.h"  /* environment_t           */
#include "k_model.h" /* k_func_t                */
#include "interp.h"  /* interp_*                */
#include "accel_k.h"

typedef struct accel_k_struct {
  interp_table_t *itable;
  k_func_t *k;
  environment_t *env;
  void *params;
} accel_k_t;

/* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
static int num_accels = 0;
static accel_k_t *accels=NULL;

/* Wrap k in the standard f(x) acceleration form */
static double k_wrap(double F, void *params)
{
  accel_k_t *p = (accel_k_t *)params;
  return (*p->k)(F, p->env, p->params);
}

static int k_tol(double FA, double kA, double FB, double kB)
{
  assert(FB > FA);
  if (FB-FA > 1e-12) {
    //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
    return 1; /* unnacceptable */
  } else {
    //printf("acceptable tol\n");
    return 0; /* acceptable */
  }
}

static int add_accel_k(k_func_t *k, environment_t *env, void *params)
{
  int i=num_accels;
  accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
  assert(accels != NULL);  
  accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
  accels[i].k = k;
  accels[i].env = env;
  accels[i].params = params;
  return i;
}

void free_accels()
{
  int i;
  for (i=0; i<num_accels; i++)
    interp_table_free(accels[i].itable);
  num_accels=0;
  if (accels) free(accels);
  accels = NULL;
}

double accel_k(k_func_t *k, double F, environment_t *env, void *params)
{
  int i;
  for (i=0; i<num_accels; i++) {
    if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
      /* We've already setup this function.
       * WARNING: we're assuming param and env are the same. */
      return interp_table_eval(accels[i].itable, accels+i, F);
    }
  }      
  /* set up a new acceleration environment */
  i = add_accel_k(k, env, params);
  return interp_table_eval(accels[i].itable, accels+i, F);
}
@ 

\section{Tension models}
\label{sec.tension_models}

TODO: tension model intro.
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.

Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
The interface is defined in [[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]].

<<tension-model.h>>=
<<license comment>>
<<tension handler types>>
<<constant tension model declarations>>
<<hooke tension model declarations>>
<<worm-like chain tension model declarations>>
<<freely-jointed chain tension model declarations>>
<<find tension definitions>>
<<tension model global declarations>>
@ 

<<tension model module makefile lines>>=
tension_model.c : sawsim.nw
        notangle -Rtension-model.c $^ > $@
tension_model.h : sawsim.nw
        notangle -Rtension-model.h $^ > $@
check_tension_model.c : sawsim.nw
        notangle -Rcheck-tension-model.c $^ > $@
check_tension_model : check_tension_model.c global.h tension_model.c tension_model.h
        gcc -g -o $@ $< tension_model.c -lcheck -lgsl -lgslcblas -lm
clean_tension_model : clean_tension_model_utils
        rm -f tension_model.c tension_model.h check_tension_model.c check_tension_model
tension_model_utils.c : sawsim.nw
        notangle -Rtension-model-utils.c $^ > $@
tension_model_utils : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
                list.c list.h tension_balance.c tension_balance.h
        gcc -g -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
tension_model_utils_static : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
                list.c list.h tension_balance.c tension_balance.h
        gcc -g -static -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
clean_tension_model_utils :
        rm -f tension_model_utils.c tension_model_utils
@ 

\subsection{Null}
\label{sec.null_tension}

For unstretchable domains.

<<null tension model getopt>>=
{NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
@ 

\subsection{Constant}
\label{sec.const_tension}

<<constant tension functions>>=
<<constant tension function>>
<<constant tension parameter create/destroy functions>>
@ 

<<constant tension model declarations>>=
<<constant tension function declaration>>
<<constant tension parameter create/destroy function declarations>>
<<constant tension model global declarations>>

@ 

An infinitely stretchable domain providing a constant tension.

<<constant tension function declaration>>=
extern double const_tension_handler(double x, void *pdata);
@ 
<<constant tension function>>=
double const_tension_handler(double x, void *pdata)
{
  tension_handler_data_t *data = (tension_handler_data_t *)pdata;
  list_t *list = data->group;
  double F;

  assert (x >= 0.0);
  list = head(list);
  assert(list != NULL); /* empty active group?! */
  F = ((const_tension_param_t *)list->d)->F;
  while (list != NULL) {
    assert(((const_tension_param_t *)list->d)->F == F);
    list = list->next;
  }
  return F;
}

@ 

<<constant tension structure>>=
typedef struct const_tension_param_struct {
  double F; /* tension (force) in N */
} const_tension_param_t;
@ 


<<constant tension parameter create/destroy function declarations>>=
extern void *string_create_const_tension_param_t(char **param_strings);
extern void destroy_const_tension_param_t(void *p);
@ 

<<constant tension parameter create/destroy functions>>=
const_tension_param_t *create_const_tension_param_t(double F)
{
  const_tension_param_t *ret
    = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
  assert(ret != NULL);
  ret->F = F;
  return ret;
}

void *string_create_const_tension_param_t(char **param_strings)
{
  return create_const_tension_param_t(atof(param_strings[0]));
}

void destroy_const_tension_param_t(void *p)
{
  if (p)
    free(p);
}

@

<<constant tension model global declarations>>=
extern int num_const_tension_params;
extern char *const_tension_param_descriptions[];
extern char const_tension_param_string[];
@ 

<<constant tension model globals>>=
int num_const_tension_params = 1;
char *const_tension_param_descriptions[] = {"tension F, N"};
char const_tension_param_string[] = "0";
@ 

<<constant tension model getopt>>=
{&const_tension_handler, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
@ 

\subsection{Hooke}
\label{sec.hooke}

<<hooke functions>>=
<<hooke function>>
<<hooke parameter create/destroy functions>>
@ 

<<hooke tension model declarations>>=
<<hooke tension function declaration>>
<<hooke tension parameter create/destroy function declarations>>
<<hooke tension model global declarations>>

@ 

The tension of a single spring is given by $F=kx$ for some spring constant $k$.
The behavior of a series of springs $k_i$ in series is given by
$$
  F = \p({\sum_i \frac{1}{k_i}})^{-1} x
$$
For a simple proof, see Appendix \ref{app.math_hooke}.

<<hooke tension function declaration>>=
extern double hooke_handler(double x, void *pdata);
@ 

<<hooke function>>=
double hooke_handler(double x, void *pdata)
{
  tension_handler_data_t *data = (tension_handler_data_t *)pdata;
  list_t *list = data->group;
  double k=0.0;

  assert (x >= 0.0);
  list = head(list);
  assert(list != NULL); /* empty active group?! */
  while (list != NULL) {
    assert( ((hooke_param_t *)list->d)->k > 0 );
    k += 1.0/ ((hooke_param_t *)list->d)->k;
    list = list->next;
  }
  k = 1.0 / k;
  return k*x;
}
@ 

<<hooke structure>>=
typedef struct hooke_param_struct {
  double k; /* spring constant in N/m */
} hooke_param_t;
@ 

<<hooke tension parameter create/destroy function declarations>>=
extern void *string_create_hooke_param_t(char **param_strings);
extern void destroy_hooke_param_t(void *p);
@ 

<<hooke parameter create/destroy functions>>=
hooke_param_t *create_hooke_param_t(double k)
{
  hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
  assert(ret != NULL);
  ret->k = k;
  return ret;
}

void *string_create_hooke_param_t(char **param_strings)
{
  return create_hooke_param_t(atof(param_strings[0]));
}

void destroy_hooke_param_t(void *p)
{
  if (p)
    free(p);
}
@

<<hooke tension model global declarations>>=
extern int num_hooke_params;
extern char *hooke_param_descriptions[];
extern char hooke_param_string[];
@ 

<<hooke tension model globals>>=
int num_hooke_params = 1;
char *hooke_param_descriptions[] = {"spring constant k, N/m"};
char hooke_param_string[]="0.05";
@ 

<<hooke tension model getopt>>=
{hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t}
@ 

\subsection{Worm-like chain}
\label{sec.wlc}

We can model several unfolded domains with a single worm-like chain.
<<worm-like chain functions>>=
<<worm-like chain function>>
<<worm-like chain handler>>
<<worm-like chain create/destroy functions>>
@ 

<<worm-like chain tension model declarations>>=
<<worm-like chain handler declaration>>
<<worm-like chain create/destroy function declarations>>
<<worm-like chain tension model global declarations>>

@ 

The combination of all unfolded domains is modeled as a worm like chain,
with the force $F_{\text{WLC}}$ approximately given by
$$
  F_{\text{WLC}} \approx \frac{k_B T}{p}
               \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
$$
where $x$ is the distance between the fixed ends,
$k_B$ is Boltzmann's constant,
$T$ is the absolute temperature,
$p$ is the persistence length, and
$L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
<<worm-like chain function>>=
double wlc(double x, double T, double p, double L)
{/* N             m         K         m         m */ 
  double xL=0.0;               /* uses global k_B */
  assert(x >= 0);
  assert(T > 0); assert(p > 0); assert(L > 0);
  if (x >= L) return HUGE_VAL;
  xL = x/L; /* unitless */
  //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
  //       k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
  return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
}
@ 
This model assumes that each unfolded domain has the same persistence length.

<<worm-like chain handler declaration>>=
extern double wlc_handler(double x, void *pdata);
@ 

<<worm-like chain handler>>=
double wlc_handler(double x, void *pdata)
{
  tension_handler_data_t *data = (tension_handler_data_t *)pdata;
  list_t *list = data->group;
  double p, L=0.0;

  list = head(list);
  assert(list != NULL); /* empty active group?! */
  p = ((wlc_param_t *)list->d)->p;
  while (list != NULL) {
    /* enforce consistent p */
    assert( ((wlc_param_t *)list->d)->p == p);
    L += ((wlc_param_t *)list->d)->L;
    list = list->next;
  }
  return wlc(x, data->env->T, p, L);
}
@ 

<<worm-like chain structure>>=
typedef struct wlc_param_struct {
  double p;      /* WLC persistence length (m) */
  double L;      /* WLC contour length (m)     */
} wlc_param_t;
@ 

<<worm-like chain create/destroy function declarations>>=
extern void *string_create_wlc_param_t(char **param_strings);
extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
@ 

<<worm-like chain create/destroy functions>>=
wlc_param_t *create_wlc_param_t(double p, double L)
{
  wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
  assert(ret != NULL);
  ret->p = p;
  ret->L = L;
  return ret;
}

void *string_create_wlc_param_t(char **param_strings)
{
  return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
}

void destroy_wlc_param_t(void *p /* wlc_param_t * */)
{
  if (p)
    free(p);
}
@

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.
TODO (now we copy the string before parsing, but still do this for clarity).
For example,
<<valid string write code>>=
char string[] = "abc";
string[1] = 'x';
@ works, but
<<valid string write code>>=
char *string = "abc";
string[1] = 'x';
@ segfaults.

<<worm-like chain tension model global declarations>>=
extern int num_wlc_params;
extern char *wlc_param_descriptions[];
extern char wlc_param_string[];
@ 
<<worm-like chain tension model globals>>=
int num_wlc_params = 2;
char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
char wlc_param_string[]="0.39e-9,28.4e-9";
@ 

<<worm-like chain tension model getopt>>=
{wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
@ 
Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.

\subsection{Freely-jointed chain}
\label{sec.fjc}

<<freely-jointed chain functions>>=
<<freely-jointed chain function>>
<<freely-jointed chain handler>>
<<freely-jointed chain create/destroy functions>>
@ 

<<freely-jointed chain tension model declarations>>=
<<freely-jointed chain handler declaration>>
<<freely-jointed chain create/destroy function declarations>>
<<freely-jointed chain tension model global declarations>>

@ 

The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
The tension of a chain of $N$ such links is given by
$$
  F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
$$
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}.
<<freely-jointed chain function>>=
double langevin(double x, void *params)
{
  if (x==0) return 0;
  return 1.0/tanh(x) - 1/x;
}
<<one dimensional bracket declaration>>
<<one dimensional solve declaration>>
double inv_langevin(double x)
{
  double lb=0.0, ub=1.0, ret;
  oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
  ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
  return ret;
}

double fjc(double x, double T, double l, int N)
{
  double L = l*(double)N;
  if (x == 0) return 0;
  //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
  return k_B*T/l * inv_langevin(x/L);
}
@ 

<<freely-jointed chain handler declaration>>=
extern double fjc_handler(double x, void *pdata);
@ 
<<freely-jointed chain handler>>=
double fjc_handler(double x, void *pdata)
{
  tension_handler_data_t *data = (tension_handler_data_t *)pdata;
  list_t *list = data->group;
  double l;
  int N=0;

  list = head(list);
  assert(list != NULL); /* empty active group?! */
  l = ((fjc_param_t *)list->d)->l;
  while (list != NULL) {
    /* enforce consistent l */
    assert( ((fjc_param_t *)list->d)->l == l);
    N += ((fjc_param_t *)list->d)->N;
    list = list->next;
  }
  return fjc(x, data->env->T, l, N);
}
@ 

<<freely-jointed chain structure>>=
typedef struct fjc_param_struct {
  double l;      /* FJC link length (m) */
  int N;         /* FJC number of links */
} fjc_param_t;
@ 

<<freely-jointed chain create/destroy function declarations>>=
extern void *string_create_fjc_param_t(char **param_strings);
extern void destroy_fjc_param_t(void *p);
@ 

<<freely-jointed chain create/destroy functions>>=
fjc_param_t *create_fjc_param_t(double l, double N)
{
  fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
  assert(ret != NULL);
  ret->l = l;
  ret->N = N;
  return ret;
}

void *string_create_fjc_param_t(char **param_strings)
{
  return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
}

void destroy_fjc_param_t(void *p)
{
  if (p)
    free(p);
}
@ 

<<freely-jointed chain tension model global declarations>>=
extern int num_fjc_params;
extern char *fjc_param_descriptions[];
extern char fjc_param_string[];
@ 

<<freely-jointed chain tension model globals>>=
int num_fjc_params = 2;
char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
char fjc_param_string[]="0.5e-9,200";
@ 

<<freely-jointed chain tension model getopt>>=
{fjc_handler, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
@ 

\subsection{Tension model implementation}

<<tension-model.c>>=
<<license comment>>
<<tension model includes>>
#include "tension_model.h"
<<tension model internal definitions>>
<<tension model globals>>
<<tension model functions>>
@ 

<<tension model includes>>=
#include <assert.h> /* assert()                */
#include <stdlib.h> /* NULL                    */
#include <math.h>   /* HUGE_VAL, sqrt(), exp() */
#include <stdio.h>  /* fprintf(), stdout       */
#include "global.h"
#include "list.h"
#include "tension_balance.h" /* oneD_*() */
@ 

<<tension model internal definitions>>=
<<constant tension structure>>
<<hooke structure>>
<<worm-like chain structure>>
<<freely-jointed chain structure>>
@ 

<<tension model globals>>=
<<tension handler array global>>
<<constant tension model globals>>
<<hooke tension model globals>>
<<worm-like chain tension model globals>>
<<freely-jointed chain tension model globals>>
@ 

<<tension model functions>>=
<<constant tension functions>>
<<hooke functions>>
<<worm-like chain functions>>
<<freely-jointed chain functions>>
@ 


\subsection{Utilities}

The tension models can get complicated, and users may want to reassure
themselves that this portion of the input physics are functioning properly. The
stand-alone program [[tension_model_utils]] demonstrates the tension model
interface without the overhead of having to understand a full simulation.
[[tension_model_utils]] takes command line model arguments like the full
simulation, and prints $F(x)$ data to the screen for the selected model over a
range of $x$.

<<tension-model-utils.c>>=
<<license comment>>
<<tension model utility includes>>
<<tension model utility definitions>>
<<tension model utility getopt functions>>
<<setup>>
<<tension model utility main>>
@ 

<<tension model utility main>>=
int main(int argc, char **argv)
{
  <<tension model getopt array definition>>
  tension_model_getopt_t *model = NULL;
  void *params;
  environment_t env;
  unsigned int flags;
  one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
  tension_handler_data_t tdata;
  int num_param_args; /* for INIT_MODEL() */
  char **param_args;  /* for INIT_MODEL() */
  get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags);
  setup(tension_handler);
  if (flags & VFLAG) {
    printf("#initializing model %s with parameters %s\n", model->name, model->params);
  }
  INIT_MODEL("utility", model, params);
  tdata.group = NULL;
  push(&tdata.group, params, 1);
  tdata.env = &env;
  tdata.persist = NULL;
  if (tension_handler[model->tg_group] == NULL) {
    printf("No tension function for model %s\n", model->name);
    exit(0);
  }
  {
    double dx=1e-10, x=0, F=0;
    printf("#F (N)\tk (%% pop. per s)\n");
    while (F >= 0 && F < 1e5 && x < 1e-6) {
      F = (*tension_handler[model->tg_group])(x, &tdata);
      printf("%g\t%g\n", x, F);
      x += dx;
    }
  }
  params = pop(&tdata.group);
  if (model->destructor)
    (*model->destructor)(params);
  return 0;
}
@ 

<<tension model utility includes>>=
#include <stdlib.h>
#include <stdio.h>
#include <string.h> /* strlen() */
#include <assert.h> /* assert() */
#include "global.h"
#include "parse.h"
#include "list.h"
#include "tension_model.h"
@ 

<<tension model utility definitions>>=
<<version definition>>
#define VFLAG 1 /* verbose */
<<string comparison definition and globals>>
<<tension model getopt definitions>>
<<initialize model definition>>

@ 

<<tension model utility getopt functions>>=
<<index tension model>>
<<tension model utility help>>
<<tension model utility get options>>
@ 

<<tension model utility help>>=
void help(char *prog_name,
          environment_t *env,
          int n_tension_models, tension_model_getopt_t *tension_models,
          int tension_model)
{
  int i, j;
  printf("usage: %s [options]\n", prog_name);
  printf("Version %s\n\n", VERSION);
  printf("A utility for understanding the available tension models\n\n");
  printf("Environmental options:\n");
  printf("-T T\tTemperature T (currently %g K)\n", env->T);
  printf("-C T\tYou can also set the temperature T in Celsius\n");
  printf("Model options:\n");
  printf("-m model\tFolded tension model (currently %s)\n",
         tension_models[tension_model].name);
  printf("-a args\tFolded tension model argument string (currently %s)\n",
         tension_models[tension_model].params);
  printf("F(x) is calculated for a range of x and printed\n");
  printf("For example:\n");
  printf("  #Distance (x)\tForce (N)\n");
  printf("  123.456\t7.89\n");
  printf("  ...\n");
  printf("-V\tChange output to verbose mode\n");
  printf("-h\tPrint this help and exit\n");
  printf("\n");
  printf("Tension models:\n");
  for (i=0; i<n_tension_models; i++) {
    printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
    for (j=0; j < tension_models[i].num_params ; j++)
      printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
    printf("\t\tdefaults: %s\n", tension_models[i].params);
  }
}
@ 

<<tension model utility get options>>=
void get_options(int argc, char **argv, environment_t *env,
                 int n_tension_models, tension_model_getopt_t *tension_models,
                 tension_model_getopt_t **model,
                 unsigned int *flags)
{
  char *prog_name = NULL;
  char c, options[] = "T:C:m:a:Vh";
  int tension_model=0, num_strings;
  extern char *optarg;
  extern int optind, optopt, opterr;

  assert (argc > 0);

  /* setup defaults */

  prog_name = argv[0];
  env->T = 300.0;   /* K           */
  *flags = 0;
  *model = tension_models;

  while ((c=getopt(argc, argv, options)) != -1) {
    switch(c) {
    case 'T':  env->T   = atof(optarg);           break;
    case 'C':  env->T   = atof(optarg)+273.15;    break;
    case 'm':
      parse_list_string(optarg, ',', NULL, NULL, &num_strings, string_array);
      tension_model = index_tension_model(n_tension_models, tension_models, optarg);
      *model = tension_models+tension_model;
      break;
    case 'a':
      tension_models[tension_model].params = optarg;
      break;
    case 'V': *flags |= VFLAG;    break;
    case '?':
      fprintf(stderr, "unrecognized option '%c'\n", optopt);
      /* fall through to default case */
    default:
      help(prog_name, env, n_tension_models, tension_models, tension_model);
      exit(1);
    }
  }
  return;
}
@ 


\section{Unfolding rate models}
\label{sec.k_models}

We have two main choices for determining $k$: Bell's model, which assumes the
folded domains exist in a single `folded' state and make instantaneous,
stochastic transitions over a free energy barrier; and Kramers' model, which
treats the unfolding as a thermalized diffusion process.
We deal with these options like we dealt with the different tension models: by bundling all model-specific 
parameters into structures, with model specific functions for determining $k$.

<<k func definition>>=
typedef double k_func_t(double F, environment_t *env, void *params);

@ 
Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
The interface is defined in [[k_model.h]], the implementation in [[k_model.c]], the unit testing in [[check_k_model.c]], and some other tidbits in [[k_model_utils.c]].

We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
because \LaTeX\ doesn't like underscores outside math mode.
Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
You could use spaces instead (`\verb+ +'),
but then your [[notangle]]s would look like [[notangle -R'k model.h']],
which I don't like as much.

<<k-model.h>>=
<<license comment>>
<<k func definition>>
<<null k declarations>>
<<const k declarations>>
<<bell k declarations>>
<<kramers k declarations>>
<<kramers integ k declarations>>
@ 

<<k model module makefile lines>>=
k_model.c : sawsim.nw
        notangle -Rk-model.c $^ > $@
k_model.h : sawsim.nw
        notangle -Rk-model.h $^ > $@
check_k_model.c : sawsim.nw
        notangle -Rcheck-k-model.c $^ > $@
check_k_model : check_k_model.c global.h k_model.c k_model.h
        gcc -g -o $@ $< k_model.c -lcheck -lgsl -lgslcblas -lm
clean_k_model : clean_k_model_utils
        rm -f k_model.c k_model.h check_k_model.c check_k_model
k_model_utils.c : sawsim.nw
        notangle -Rk-model-utils.c $^ > $@
k_model_utils : k_model_utils.c global.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h
        gcc -g -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
k_model_utils_static : k_model_utils.c global.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h
        gcc -g -static -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
clean_k_model_utils :
        rm -f k_model_utils.c k_model_utils
@ 

\subsection{Null}
\label{sec.null_k}

For domains that are never folded.

<<null k declarations>>=
@ 
<<null k functions>>=
@ 
<<null k globals>>=
@ 

<<null k model getopt>>=
{"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
@ 

\subsection{Constant rate model}
\label{sec.k_const}

This is a pretty boring model, but it is usefull for testing the engine.
$$
  k = k_0\;,
$$
where $k_0$ is the constant unfolding rate.

<<const k functions>>=
<<const k function>>
<<const k structure create/destroy functions>>
@ 

<<const k declarations>>=
<<const k function declaration>>
<<const k structure create/destroy function declarations>>
<<const k global declarations>>

@ 

<<const k structure>>=
typedef struct const_k_param_struct {
  double knot;   /* transition rate k_0 (frac population per s) */
} const_k_param_t;
@ 

<<const k function declaration>>=
double const_k(double F, environment_t *env, void *const_k_params);
@ 
<<const k function>>=
double const_k(double F, environment_t *env, void *const_k_params)
{ /* Returns the rate constant k in frac pop per s. */
  const_k_param_t *p = (const_k_param_t *) const_k_params;
  assert(p != NULL);
  assert(p->knot > 0);
  return p->knot;
}
@ 

<<const k structure create/destroy function declarations>>=
void *string_create_const_k_param_t(char **param_strings);
void destroy_const_k_param_t(void *p);
@ 

<<const k structure create/destroy functions>>=
const_k_param_t *create_const_k_param_t(double knot)
{
  const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
  assert(ret != NULL);
  ret->knot = knot;
  return ret;
}

void *string_create_const_k_param_t(char **param_strings)
{
  return create_const_k_param_t(atof(param_strings[0]));
}

void destroy_const_k_param_t(void *p /* const_k_param_t * */)
{
  if (p)
    free(p);
}
@

<<const k global declarations>>=
extern int num_const_k_params;
extern char *const_k_param_descriptions[];
extern char const_k_param_string[];
@

<<const k globals>>=
int num_const_k_params = 1;
char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
char const_k_param_string[]="1";
@

<<const k model getopt>>=
{"const", "constant rate transitions", &const_k, num_const_k_params, const_k_param_descriptions, (char *)const_k_param_string, &string_create_const_k_param_t, &destroy_const_k_param_t}
@ 

\subsection{Bell's model}
\label{sec.bell}

Quantitatively, Bell's model gives $k$ as 
$$
  k = k_0 \cdot e^{F dx / k_B T} \;,
$$
where $k_0$ is the unforced unfolding rate,
$F$ is the applied unfolding force,
$dx$ is the distance to the transition state, and
$k_B T$ is the thermal energy\citep{schlierf06}.

<<bell k functions>>=
<<bell k function>>
<<bell k structure create/destroy functions>>
@ 

<<bell k declarations>>=
<<bell k function declaration>>
<<bell k structure create/destroy function declarations>>
<<bell k global declarations>>

@ 

<<bell k structure>>=
typedef struct bell_param_struct {
  double knot;   /* transition rate k_0 (frac population per s) */
  double dx;     /* `distance to transition state' (nm) */
} bell_param_t;
@ 

<<bell k function declaration>>=
double bell_k(double F, environment_t *env, void *bell_params);
@ 
<<bell k function>>=
double bell_k(double F, environment_t *env, void *bell_params)
{ /* Returns the rate constant k in frac pop per s.
   * Takes F in N, T in K, knot in frac pop per s, and dx in m.
   * uses global k_B in J/K */
  bell_param_t *p = (bell_param_t *) bell_params;
  assert(F >= 0); assert(env->T > 0);
  assert(p != NULL);
  assert(p->knot > 0); assert(p->dx >= 0);
/*
  printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
         F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
         p->knot * exp(F*p->dx / (k_B*env->T) ));
  printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
*/
  return p->knot * exp(F*p->dx / (k_B*env->T) );
}
@ 

<<bell k structure create/destroy function declarations>>=
void *string_create_bell_param_t(char **param_strings);
void destroy_bell_param_t(void *p);
@ 

<<bell k structure create/destroy functions>>=
bell_param_t *create_bell_param_t(double knot, double dx)
{
  bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
  assert(ret != NULL);
  ret->knot = knot;
  ret->dx = dx;
  return ret;
}

void *string_create_bell_param_t(char **param_strings)
{
  return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
}

void destroy_bell_param_t(void *p /* bell_param_t * */)
{
  if (p)
    free(p);
}
@

<<bell k global declarations>>=
extern int num_bell_params;
extern char *bell_param_descriptions[];
extern char bell_param_string[];
@

<<bell k globals>>=
int num_bell_params = 2;
char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
char bell_param_string[]="3.3e-4,0.25e-9";
@

<<bell k model getopt>>=
{"bell", "thermalized, two-state transitions", &bell_k, num_bell_params, bell_param_descriptions, (char *)bell_param_string, &string_create_bell_param_t, &destroy_bell_param_t}
@ 
Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696

\subsection{Kramer's model}
\label{sec.kramers}

Kramer's model gives $k$ as
%$$
%  \frac{1}{k} = \frac{1}{D}
%                \int_{x_\text{min}}^{x_\text{max}}
%                     e^{\frac{-U_F(x)}{k_B T}}
%                     \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
%                     dx,
%$$
%where $D$ is the diffusion constant,
%$U_F(x)$ is the free energy along the unfolding coordinate, and
%$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
% before the minimum and after the maximum, respectively \citep{schlierf06}.
$$
  k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
$$
where $D$ is the diffusion constant,
$$
  l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
$$
are length scales where
$x_c(F)$ is the location of the bound state and
$x_{ts}(F)$ is the location of the transition state,
$E(F,x)$ is the free energy, and
$E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
(\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
 \citet{evans97} Eqn.~3).

<<kramers k functions>>=
<<kramers k function>>
<<kramers k structure create/destroy functions>>
@ 

<<kramers k declarations>>=
<<kramers k function declaration>>
<<kramers k structure create/destroy function declarations>>
<<kramers k global declarations>>

@ 

<<kramers k structure>>=
//typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
//typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);

typedef struct kramers_param_struct {
  double D;                 /* diffusion rate (um/s)                 */
  char *xc;
  char *xts;
  char *ddEdxx;
  char *E;
  FILE *in;
  FILE *out;
  //kramers_x_func_t *xc;     /* function returning metastable x (nm)  */
  //kramers_x_func_t *xts;    /* function returning transition x (nm)  */
  //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
  //kramers_E_func_t *E;      /* function returning E (J)              */
  //double *E_params;         /* parametrize the energy landscape     */
  //int n_E_params;           /* length of E_params array              */
} kramers_param_t;
@ 

<<kramers k function declaration>>=
extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
extern double kramers_k(double F, environment_t *env, void *kramers_params);
@ 
<<kramers k function>>=
double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
{
  static char input[160]={0};
  static char output[80]={0};
  /* setup the environment */
  fprintf(in, "F = %g; T = %g\n", F, T);
  sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
  string_eval(in, out, input, 80, output);
  return atof(output);
}

double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
{
  static char input[160]={0};
  static char output[80]={0};
  /* setup the environment */
  fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
  sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
  string_eval(in, out, input, 80, output);
  return atof(output);
}

double kramers_E(double F, environment_t *env, void *kramers_params, double x)
{
  kramers_param_t *p = (kramers_param_t *) kramers_params;
  return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
}

double kramers_k(double F, environment_t *env, void *kramers_params)
{ /* Returns the rate constant k in frac pop per s.
   * Takes F in N, T in K, knot in frac pop per s, and dx in m.
   * uses global k_B in J/K */
  kramers_param_t *p = (kramers_param_t *) kramers_params;
  double kbT, xc, xts, lc, lts, Eb;
  assert(F >= 0); assert(env->T > 0);
  assert(p != NULL);
  assert(p->D > 0);
  assert(p->xc != NULL); assert(p->xts != NULL);
  assert(p->ddEdxx != NULL); assert(p->E != NULL);
  assert(p->in != NULL); assert(p->out != NULL);
  kbT = k_B*env->T;
  xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
  if (xc == -1.0) return -1.0; /* error (Too much force) */
  xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
  if (xc == -1.0) return -1.0; /* error (Too much force) */
  lc = sqrt(2.0*M_PI*kbT /
            kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
  lts = sqrt(-2.0*M_PI*kbT/
            kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
  Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
     - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
  //fprintf(stderr,"D = %g, lc = %g, lts = %g, Eb = %g, kbT = %g, k = %g\n", p->D, lc, lts, Eb, kbT, p->D/(lc*lts) * exp(-Eb/kbT));
  return p->D/(lc*lts) * exp(-Eb/kbT);
}
@ 

<<kramers k structure create/destroy function declarations>>=
void *string_create_kramers_param_t(char **param_strings);
void destroy_kramers_param_t(void *p);
@ 

<<kramers k structure create/destroy functions>>=
kramers_param_t *create_kramers_param_t(double D,
                                        char *xc_fn, char *xts_fn,
                                        char *E_fn, char *ddEdxx_fn,
                                        char *global_define_string)
//                              kramers_x_func_t *xc, kramers_x_func_t *xts,
//                            kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
//                            double *E_params, int n_E_params)
{
  kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
  assert(ret != NULL);
  ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
  assert(ret->xc != NULL);
  ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
  assert(ret->xts != NULL);
  ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
  assert(ret->E != NULL);
  ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
  assert(ret->ddEdxx != NULL);
  ret->D = D;
  strcpy(ret->xc, xc_fn);
  strcpy(ret->xts, xts_fn);
  strcpy(ret->E, E_fn);
  strcpy(ret->ddEdxx, ddEdxx_fn);
  //ret->E_params = E_params;
  //ret->n_E_params = n_E_params;
  string_eval_setup(&ret->in, &ret->out);
  fprintf(ret->in, "kB = %g\n", k_B);
  fprintf(ret->in, "%s\n", global_define_string);
  return ret;
}

/* takes an array of 4 functions: {"(1,2),(2,0),..."} */
void *string_create_kramers_param_t(char **param_strings)
{
  return create_kramers_param_t(atof(param_strings[0]),
                                param_strings[2],
                                param_strings[3],
                                param_strings[4],
                                param_strings[5],
                                param_strings[1]);
}

void destroy_kramers_param_t(void *p /* kramers_param_t * */)
{
  kramers_param_t *pk = (kramers_param_t *)p;
  if (pk) {
    string_eval_teardown(&pk->in, &pk->out);
    //if (pk->E_params)
    //  free(pk->E_params);
    free(pk);
  }
}
@

For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
Schlierf and Rief used a cubic-spline interpolation routine and the double integral form of Kramers' theory, so we get to pick an actual function to approximate the energy landscape.
We follow \citet{shillcock98} and use
$$ 
  E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
                         \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
                        -F x
$$
where TODO, variable meanings.%\citep{shillcock98}.
The first and second derivatives of $E(F,E_b,x,x_s)$ are
\begin{align}
  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\\
  E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
\end{align}
$E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
$b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
The roots are therefor at
\begin{align}
  x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
        &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
        &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
\end{align}

<<kramers k global declarations>>=
extern int num_kramers_params;
extern char *kramers_param_descriptions[];
extern char kramers_param_string[];
@

<<kramers k globals>>=
int num_kramers_params = 6;
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"};
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)}";
@
Which we initialized with parameters for ddFLN4\citep{schlierf06}.
There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
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.
The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
It works so long as [[val_1]] is not [[false]].

<<kramers k model getopt>>=
{"kramers", "thermalized, diffusion-limited transitions (function parameters are parsed by Python.  For distances, the environment variables force 'F' and temperature 'T' are defined, as well as the Boltzmann constant 'kB'.  For energies, the position 'x' is also defined.  Functions must evaluate to a float representing their output, and produce no output on stdout.", &kramers_k, num_kramers_params, kramers_param_descriptions, (char *)kramers_param_string, &string_create_kramers_param_t, &destroy_kramers_param_t}
@ 

\subsection{Kramer's model (natural cubic splines)}
\label{sec.kramers_integ}

Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
(\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
\citet{schlierf06} Eqn.~2)
$$
  \frac{1}{k} = \frac{1}{D}
                \int_{x_\text{min}}^{x_\text{max}}
                     e^{\frac{U_F(x)}{k_B T}}
                     \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
                     dx,
$$
where $D$ is the diffusion constant,
$U_F(x)$ is the free energy along the unfolding coordinate, and
$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.

We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
interpolating between them with cubic splines.
We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.

<<kramers integ k functions>>=
<<spline functions>>
<<kramers integ A k functions>>
<<kramers integ B k functions>>
<<kramers integ k function>>
<<kramers integ k structure create/destroy functions>>
@ 

<<kramers integ k declarations>>=
<<kramers integ k function declaration>>
<<kramers integ k structure create/destroy function declarations>>
<<kramers integ k global declarations>>

@ 

<<kramers integ k structure>>=
typedef double func_t(double x, void *params);
typedef struct kramers_integ_param_struct {
  double D;            /* diffusion rate (um/s)                 */
  double x_min;        /* integration bounds                    */
  double x_max;
  func_t *E;           /* function returning E (J)              */
  void *E_params;      /* parametrize the energy landscape     */
  destroy_data_func_t *destroy_E_params;

  double F;            /* for passing into GSL evaluations      */
  environment_t *env;
} kramers_integ_param_t;
@ 

<<spline param structure>>=
typedef struct spline_param_struct {
  int n_knots;            /* natural cubic spline knots for U(x)   */
  double *x;
  double *y;
  gsl_spline *spline;    /* GSL spline parameters               */
  gsl_interp_accel *acc; /* GSL spline acceleration data        */
} spline_param_t;
@ 

We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
$$
  \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
$$
<<kramers integ A k functions>>=
double kramers_integ_k_integrandA(double x, void *params)
{
  kramers_integ_param_t *p = (kramers_integ_param_t *)params;
  double U = (*p->E)(x, p->E_params);
  double Fx = p->F * x;
  double kbT = k_B * p->env->T;
  //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
  return exp(-(U-Fx)/kbT);
}
double kramers_integ_k_integralA(double x, void *params)
{
  kramers_integ_param_t *p = (kramers_integ_param_t *)params;
  gsl_function f;
  double result, abserr;
  size_t neval; /* for qng */
  static gsl_integration_workspace *w = NULL;
  if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
  f.function = &kramers_integ_k_integrandA;
  f.params = params;
  //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
  assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
  //fprintf(stderr, "integralA = %g\n", result);
  return result;
}
@ 

$$
  \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
                \int_{x_\text{min}}^{x_\text{max}}
                     e^{\frac{U_F(x)}{k_B T}}
                     \text{Integral}_A(x)
                     dx,
$$
<<kramers integ B k functions>>=
double kramers_integ_k_integrandB(double x, void *params)
{
  kramers_integ_param_t *p = (kramers_integ_param_t *)params;
  double U = (*p->E)(x, p->E_params);
  double Fx = p->F * x;
  double kbT = k_B * p->env->T;
  double intA = kramers_integ_k_integralA(x, params);
  //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
  return exp((U-Fx)/kbT)*intA;
}
double kramers_integ_k_integralB(double F, environment_t *env, void *params)
{
  kramers_integ_param_t *p = (kramers_integ_param_t *)params;
  gsl_function f;
  double result, abserr;
  size_t neval; /* for qng */
  static gsl_integration_workspace *w = NULL;
  if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
  f.function = &kramers_integ_k_integrandB;
  f.params = params;
  p->F = F;
  p->env = env;
  //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
  assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
  //fprintf(stderr, "integralB = %g\n", result);
  return result;
}
@ 

With the integrals taken care of, evaluating $k$ is simply
$$
  k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
$$
<<kramers integ k function declaration>>=
double kramers_integ_k(double F, environment_t *env, void *kramers_params);
@ 
<<kramers integ k function>>=
double kramers_integ_k(double F, environment_t *env, void *kramers_params)
{ /* Returns the rate constant k in frac pop per s. Takes F in N. */
  kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
  return p->D/kramers_integ_k_integralB(F, env, kramers_params);
}
@ 

<<kramers integ k structure create/destroy function declarations>>=
void *string_create_kramers_integ_param_t(char **param_strings);
void destroy_kramers_integ_param_t(void *p);
@ 

<<kramers integ k structure create/destroy functions>>=
kramers_integ_param_t *create_kramers_integ_param_t(double D,
                              double xmin, double xmax, func_t *E,
                              void *E_params,
                              destroy_data_func_t *destroy_E_params)
{
  kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
  assert(ret != NULL);
  ret->D = D;
  ret->x_min = xmin;
  ret->x_max = xmax;
  ret->E = E;
  ret->E_params = E_params;
  ret->destroy_E_params = destroy_E_params;
  return ret;
}

void *string_create_kramers_integ_param_t(char **param_strings)
{
  //printf("create kramers integ.  D: %s, knots: %s, x in [%s, %s]\n",
  //       param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
  void *E_params = string_create_spline_param_t(param_strings+1);
  return create_kramers_integ_param_t(atof(param_strings[0]),
                                      atof(param_strings[2]),
                                      atof(param_strings[3]),
                                      &spline_eval, E_params,
                                      destroy_spline_param_t);
}

void destroy_kramers_integ_param_t(void *params)
{
  kramers_integ_param_t *p = (kramers_integ_param_t *)params;
  if (p) {
    if (p->E_params)
      (*p->destroy_E_params)(p->E_params);
    free(p);
  }
}
@

Finally we have the GSL spline wrappers:
<<spline functions>>=
<<spline function>>
<<create/destroy spline>>
@ 

<<spline function>>=
double spline_eval(double x, void *spline_params)
{
  spline_param_t *p = (spline_param_t *)(spline_params);
  return gsl_spline_eval(p->spline, x, p->acc);
}
@ 

<<create/destroy spline>>=
spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
{
  spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
  assert(ret != NULL);
  ret->n_knots = n_knots;
  ret->x = x;
  ret->y = y;
  ret->acc = gsl_interp_accel_alloc();
  ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
  gsl_spline_init(ret->spline, x, y, n_knots);
  return ret;
}

/* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
void *string_create_spline_param_t(char **param_strings)
{
  char **knot_coords;
  int i, num_knots;
  double *x=NULL, *y=NULL;
  /* break into ordered pairs */
  parse_list_string(param_strings[0], SEP, '(', ')',
                    &num_knots, &knot_coords);
  x = (double *)malloc(sizeof(double)*num_knots);
  assert(x != NULL);
  y = (double *)malloc(sizeof(double)*num_knots);
  assert(x != NULL);
  for (i=0; i<num_knots; i++) {
    if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
      fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
      exit(1);
    }
    //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
  }
  free_parsed_list(num_knots, knot_coords);
  return create_spline_param_t(num_knots, x, y);
}

void destroy_spline_param_t(void *params)
{
  spline_param_t *p = (spline_param_t *)params;
  if (p) {
    if (p->spline)
      gsl_spline_free(p->spline);
    if (p->acc)
      gsl_interp_accel_free(p->acc);
    if (p->x)
      free(p->x);
    if (p->y)
      free(p->y);
    free(p);
  }
}
@ 

<<kramers integ k global declarations>>=
extern int num_kramers_integ_params;
extern char *kramers_integ_param_descriptions[];
extern char kramers_integ_param_string[];
@

<<kramers integ k globals>>=
int num_kramers_integ_params = 4;
char *kramers_integ_param_descriptions[] = {
                           "diffusion D, m^2 / s",
                           "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
                           "starting integration bound x_min, m",
                           "ending integration bound x_max, m"};
//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";
//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";
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";
/* D from Schlierf 2006, Biophys 90, 4, L33, sup.
 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
@

<<kramers integ k model getopt>>=
{"kramers_integ", "thermalized, diffusion-limited transitions", &kramers_integ_k, num_kramers_integ_params, kramers_integ_param_descriptions, (char *)kramers_integ_param_string, &string_create_kramers_integ_param_t, &destroy_kramers_integ_param_t}
@ 
Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696

\subsection{Unfolding model implementation}

<<k-model.c>>=
<<license comment>>
<<k model includes>>
#include "k_model.h"
<<k model internal definitions>>
<<k model globals>>
<<k model functions>>
@ 

<<k model includes>>=
#include <assert.h> /* assert()                */
#include <stdlib.h> /* NULL, malloc()          */
#include <math.h>   /* HUGE_VAL, sqrt(), exp() */
#include <stdio.h>  /* fprintf(), stdout       */
#include <string.h> /* strlen(), strcpy()      */
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include "global.h"
#include "parse.h"
@ 

<<k model internal definitions>>=
<<const k structure>>
<<bell k structure>>
<<kramers k structure>>
<<kramers integ k structure>>
<<spline param structure>>
@ 

<<k model globals>>=
<<null k globals>>
<<const k globals>>
<<bell k globals>>
<<kramers k globals>>
<<kramers integ k globals>>
@ 

<<k model functions>>=
<<null k functions>>
<<const k functions>>
<<bell k functions>>
<<kramers k functions>>
<<kramers integ k functions>>
@ 

\subsection{Unfolding model unit tests}

Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-k-model.c>>=
<<k model check includes>>
<<check relative error>>
<<model globals>>
<<k model test suite>>
<<main check program>>
@ 

<<k model check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
#include <stdio.h>  /* sprintf()                             */
#include <assert.h> /* assert()                              */
#include <math.h>   /* exp()                                 */
<<check includes>>
#include "global.h"
#include "k_model.h"
@ 

<<k model test suite>>=
<<const k tests>>
<<bell k tests>>
<<k model suite definition>>
@ 

<<k model suite definition>>=
Suite *test_suite (void)
{
  Suite *s = suite_create ("k model");
  <<const k test case defs>>
  <<bell k test case defs>>

  <<const k test case adds>>
  <<bell k test case adds>>
  return s;
}
@ 

\subsubsection{Constant}

<<const k test case defs>>=
TCase *tc_const_k = tcase_create("Constant k");
@ 

<<const k test case adds>>=
tcase_add_test(tc_const_k, test_const_k_create_destroy);
tcase_add_test(tc_const_k, test_const_k_over_range);
suite_add_tcase(s, tc_const_k);
@ 


<<const k tests>>=
START_TEST(test_const_k_create_destroy)
{
  char *knot[] = {"1","2","3","4","5","6"};
  char *params[] = {knot[0]};
  void *p = NULL;
  int i;
  for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
    params[0] = knot[i];
    p = string_create_const_param_t(params); 
    destroy_const_param_t(p);
  }
}
END_TEST

START_TEST(test_const_k_over_range)
{
  double F, Fm=0.0, FM=10, dF=.1, T=300.0;
  char *knot[] = {"1","2","3","4","5","6"};
  char *params[] = {knot[0]};
  void *p = NULL;
  environment_t env;
  char param_string[80];
  int i;
  env.T = T;
  for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
    params[0] = knot[i];
    p = string_create_const_param_t(params); 
    for ( F=Fm, F<FM, F+=dF ) {
      fail_unless(const_k(F, &env, p)==atof(knot[i]),
          "const_k(%g, %g, &{%s,%s}) = %g != %s",
          F, T, knot[i], dx, const_k(F, &env, p), knot[i]);
    }
    destroy_const_param_t(p);
  }
}
END_TEST
@ 

\subsubsection{Bell}

<<bell k test case defs>>=
TCase *tc_bell = tcase_create("Bell's k");
@ 

<<bell k test case adds>>=
tcase_add_test(tc_bell, test_bell_k_create_destroy);
tcase_add_test(tc_bell, test_bell_k_at_zero);
tcase_add_test(tc_bell, test_bell_k_at_one);
suite_add_tcase(s, tc_bell);
@ 

<<bell k tests>>=
START_TEST(test_bell_k_create_destroy)
{
  char *knot[] = {"1","2","3","4","5","6"};
  char *dx="1";
  char *params[] = {knot[0], dx};
  void *p = NULL;
  int i;
  for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
    params[0] = knot[i];
    p = string_create_bell_param_t(params); 
    destroy_bell_param_t(p);
  }
}
END_TEST

START_TEST(test_bell_k_at_zero)
{
  double F=0.0, T=300.0;
  char *knot[] = {"1","2","3","4","5","6"};
  char *dx="1";
  char *params[] = {knot[0], dx};
  void *p = NULL;
  environment_t env;
  char param_string[80];
  int i;
  env.T = T;
  for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
    params[0] = knot[i];
    p = string_create_bell_param_t(params); 
    fail_unless(bell_k(F, &env, p)==atof(knot[i]),
                "bell_k(%g, %g, &{%s,%s}) = %g != %s",
                F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
    destroy_bell_param_t(p);
  }
}
END_TEST

START_TEST(test_bell_k_at_one)
{
  double T=300.0;
  char *knot[] = {"1","2","3","4","5","6"};
  char *dx="1";
  char *params[] = {knot[0], dx};
  double F= k_B*T/atof(dx);
  void *p = NULL;
  environment_t env;
  char param_string[80];
  int i;
  env.T = T;
  for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
    params[0] = knot[i];
    p = string_create_bell_param_t(params); 
    CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
    //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
    //            "bell_k(%g, %g, &{%s,%s}) = %g != %g",
    //            F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
    destroy_bell_param_t(p);
  }
}
END_TEST
@ 

<<kramers k tests>>=
@ 

<<kramers k test case def>>=
@ 

<<kramers k test case add>>=
@ 

<<k model function tests>>=
<<check relative error>>

START_TEST(test_something)
{
  double k=5, x=3, last_x=2.0, F;
  one_dim_fn_t *handlers[] = {&hooke};
  void *data[] =               {&k};
  double xi[] =                {0.0};
  int active[] =               {1};
  int new_active_groups = 1;
  F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
  fail_unless(F = k*x, NULL);
}
END_TEST
@ 

\subsection{Utilities}

The unfolding models can get complicated, and are sometimes hard to explain to others.
The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
the overhead of having to understand a full simulation.
[[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
or special mode, where the operation depends on the specific model selected.

<<k-model-utils.c>>=
<<license comment>>
<<k model utility includes>>
<<k model utility definitions>>
<<k model utility getopt functions>>
<<k model utility multi model E>>
<<k model utility main>>
@ 

<<k model utility main>>=
int main(int argc, char **argv)
{
  <<k model getopt array definition>>
  k_model_getopt_t *model = NULL;
  void *params;
  enum mode_t mode;
  environment_t env;
  unsigned int flags;
  int num_param_args; /* for INIT_MODEL() */
  char **param_args;  /* for INIT_MODEL() */
  double Fmax;
  double special_xmin;
  double special_xmax;
  get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
              &Fmax, &special_xmin, &special_xmax, &flags);
  if (flags & VFLAG) {
    printf("#initializing model %s with parameters %s\n", model->name, model->params);
  }
  INIT_MODEL("utility", model, params);
  switch (mode) {
    case M_K_OF_F :
      if (model->k == NULL) {
        printf("No k function for model %s\n", model->name);
        exit(0);
      }
      if (Fmax <= 0) {
        printf("Fmax = %g <= 0\n", Fmax);
        exit(1);
      }
      {
        double F, k=1.0;
        int i, N=200;
        printf("#F (N)\tk (%% pop. per s)\n");
        for (i=0; i<=N; i++) {
          F = Fmax*i/(double)N;
          k = (*model->k)(F, &env, params);
          if (k < 0) break;
          printf("%g\t%g\n", F, k);
        }
      }
      break;
    case M_SPECIAL :
      if (model->k == NULL || model->k == &bell_k) {
        printf("No special mode for model %s\n", model->name);
        exit(1);
      }
      if (special_xmax <= special_xmin) {
        printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
        exit(1);
      }
      {
        double Xrng=(special_xmax-special_xmin), x, E;
        int i, N=500;
        printf("#x (m)\tE (J)\n");
        for (i=0; i<=N; i++) {
          x = special_xmin + Xrng*i/(double)N;
          E = multi_model_E(model->k, params, &env, x);
          printf("%g\t%g\n", x, E);
        }
      }
      break;
    default :
      break;
  }
  if (model->destructor)
    (*model->destructor)(params);
  return 0;
}
@ 

<<k model utility multi model E>>=
double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
{
  kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
  double E;
  if (k == kramers_integ_k)
    E = (*p->E)(x, p->E_params);
  else 
    E = kramers_E(0, env, model_params, x);
  return E;
}

@ 

<<k model utility includes>>=
#include <stdlib.h>
#include <stdio.h>
#include <string.h> /* strlen() */
#include <assert.h> /* assert() */
#include "global.h"
#include "parse.h"
#include "string_eval.h"
#include "k_model.h"
@ 

<<k model utility definitions>>=
<<version definition>>
#define VFLAG 1 /* verbose */
enum mode_t {M_K_OF_F, M_SPECIAL};
<<string comparison definition and globals>>
<<k model getopt definitions>>
<<initialize model definition>>
<<kramers k structure>>
<<kramers integ k structure>>

@ 

<<k model utility getopt functions>>=
<<index k model>>
<<k model utility help>>
<<k model utility get options>>
@ 

<<k model utility help>>=
void help(char *prog_name,
          environment_t *env,
          int n_k_models, k_model_getopt_t *k_models,
          int k_model)
{
  int i, j;
  printf("usage: %s [options]\n", prog_name);
  printf("Version %s\n\n", VERSION);
  printf("A utility for understanding the available unfolding models\n\n");
  printf("Environmental options:\n");
  printf("-T T\tTemperature T (currently %g K)\n", env->T);
  printf("-C T\tYou can also set the temperature T in Celsius\n");
  printf("Model options:\n");
  printf("-k model\tTransition rate model (currently %s)\n",
         k_models[k_model].name);
  printf("-K args\tTransition rate model argument string (currently %s)\n",
         k_models[k_model].params);
  printf("There are two output modes.  In standard mode, k(F) is printed\n");
  printf("For example:\n");
  printf("  #Force (N)\tk (% pop. per s)\n");
  printf("  123.456\t7.89\n");
  printf("  ...\n");
  printf("In special mode, the output depends on the model.\n");
  printf("For models defining an energy function E(x), that function is printed\n");
  printf("  #Position (m)\tE (J)\n");
  printf("  0.0012\t0.0034\n");
  printf("  ...\n");
  printf("-m\tChange output to standard mode\n");
  printf("-M\tChange output to special mode\n");
  printf("-F\tSet the maximum F value for the standard mode k(F) output\n");
  printf("-x\tSet the minimum x value for the special mode E(x) output\n");
  printf("-X\tSet the maximum x value for the special mode E(x) output\n");
  printf("-V\tChange output to verbose mode\n");
  printf("-h\tPrint this help and exit\n");
  printf("\n");
  printf("Unfolding rate models:\n");
  for (i=0; i<n_k_models; i++) {
    printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
    for (j=0; j < k_models[i].num_params ; j++)
      printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
    printf("\t\tdefaults: %s\n", k_models[i].params);
  }
}
@ 

<<k model utility get options>>=
void get_options(int argc, char **argv, environment_t *env,
                 int n_k_models, k_model_getopt_t *k_models,
                 enum mode_t *mode, k_model_getopt_t **model,
                 double *Fmax, double *special_xmin, double *special_xmax,
                 unsigned int *flags)
{
  char *prog_name = NULL;
  char c, options[] = "T:C:k:K:mMF:x:X:Vh";
  int k_model=0;
  extern char *optarg;
  extern int optind, optopt, opterr;

  assert (argc > 0);

  /* setup defaults */

  prog_name = argv[0];
  env->T = 300.0;   /* K           */
  *mode = M_K_OF_F;
  *flags = 0;
  *model = k_models;
  *Fmax = 1e-9;
  *special_xmax = 1e-8;
  *special_xmin = 0.0;

  while ((c=getopt(argc, argv, options)) != -1) {
    switch(c) {
    case 'T':  env->T   = atof(optarg);           break;
    case 'C':  env->T   = atof(optarg)+273.15;    break;
    case 'k':
      k_model = index_k_model(n_k_models, k_models, optarg);
      *model = k_models+k_model;
      break;
    case 'K':
      k_models[k_model].params = optarg;
      break;
    case 'm': *mode = M_K_OF_F;             break;
    case 'M': *mode = M_SPECIAL;            break;
    case 'F': *Fmax = atof(optarg);         break;
    case 'x': *special_xmin = atof(optarg); break;
    case 'X': *special_xmax = atof(optarg); break;
    case 'V': *flags |= VFLAG;              break;
    case '?':
      fprintf(stderr, "unrecognized option '%c'\n", optopt);
      /* fall through to default case */
    default:
      help(prog_name, env, n_k_models, k_models, k_model);
      exit(1);
    }
  }
  return;
}
@ 


\section{\LaTeX\ documentation}

The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
The comment blocks are extracted (with nicely formatted code blocks), using
<<latex makefile lines>>=
sawsim.tex : sawsim.nw
        noweave -latex -index -delay $^ > $@
@ 
and compiled using
<<latex makefile lines>>=
sawsim.pdf : sawsim.tex sawsim.bib
        pdflatex $^
        bibtex sawsim
        pdflatex $^
        pdflatex $^
@
The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.

[[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
<<latex makefile lines>>=
semi-clean_latex : 
        rm -f sawsim.aux sawsim.bbl sawsim.blg sawsim.log sawsim.out sawsim.tex
clean_latex : semi-clean_latex
        rm -f sawsim.pdf
@ 

\section{Makefile}

[[make]] is a common utility on *nix systems for managing dependencies while building software.
In this case, we have several \emph{targets} that we'd like to build:
 the main simulation program \prog;
 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
 the documentation [[sawsim.pdf]];
 and the [[Makefile]] itself.
Besides the generated files, there is the utility target
 [[clean]] for removing all generated files except [[Makefile]].
% and
% [[dist]] for generating a distributable tar file.

Extract the makefile with
`[[notangle -Rmakefile sawsim.nw | sed 's/        /\t/' > Makefile]]'.
The sed is needed because notangle replaces tabs with eight spaces,
but [[make]] needs tabs.
<<makefile>>=
all : sawsim sawsim.pdf Makefile tension_model_utils k_model_utils

view : sawsim.pdf
        xpdf $^ &
profile : sawsim_profile
        sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ -mnull -Mwlc -A0.39e-9,28e-9 -F8        
        gprof sawsim_profile gmon.out > $@

clean : clean_check clean_list clean_tension clean_k_model clean_parse clean_string_eval clean_latex
        rm -f sawsim.c sawsim sawsim_static sawsim_profile gmon.out global.h

sawsim.c : sawsim.nw
        notangle $^ > $@
global.h : sawsim.nw
        notangle -Rglobal.h $^ > $@
sawsim : sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
                tension_model.c tension_model.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h \
                accel_k.c accel_k.h
        gcc -g -o $@ $< list.c tension_balance.c tension_model.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm
sawsim_static : sawsim
        gcc -g -static -o $@ sawsim.c list.c tension_balance.c tension_model.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm
sawsim_profile : sawsim
        gcc -g -pg -o $@ sawsim.c list.c tension_balance.c tension_model.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm


check_sawsim.c : sawsim.nw
        notangle -Rchecks $^ > $@
check_sawsim : check_sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
                k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h accel_k.c accel_k.h
        gcc -g -o $@ $< list.c tension_balance.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm -lcheck
clean_check :
        rm -f check_sawsim check_sawsim.c
check : check_list check_tension_balance check_k_model check_parse check_string_eval check_accel_k check_sawsim
        check_list
        check_tension_balance
        check_k_model
        check_parse
        check_string_eval
        check_accel_k
        check_sawsim

<<list module makefile lines>>
<<tension balance module makefile lines>>
<<tension model module makefile lines>>
<<k model module makefile lines>>
<<parse module makefile lines>>
<<string eval module makefile lines>>
<<accel k module makefile lines>>
<<latex makefile lines>>

Makefile : sawsim.nw
        notangle -Rmakefile $^ | sed 's/        /\t/' > $@
@ 
This is fairly self-explanatory.
We compile a dynamically linked version ([[sawsim]]) and a statically linked version ([[sawsim_static]]).
The static version is about 9 times as big, but it runs without needing dynamic GSL libraries to link to, so it's a better format for distributing to a cluster for parallel evaluation.


\section{Math}

\subsection{Hookean springs in series}
\label{app.math_hooke}

\begin{theorem}
The effective spring constant for $n$ springs $k_i$ in series is given by
$$
  k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
$$
\end{theorem}

\begin{proof}
\begin{align}
  F &= k_i x_i &&\forall i \le n \\
  x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
  x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
  F &= k_1 x_1 = k_\text{eff} x \\
  k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}} 
                = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
\end{align}
\end{proof}

\phantomsection
\addcontentsline{toc}{section}{References}
\bibliography{sawsim}

\end{document}