testing/const_rate/const_rate.sh

   1 #!/bin/bash
   2 #
   3 # Use a single hooke domain with a constant unfolding rate, so unfolding force
   4 # is proportional to time, and we expect an exponential population decay.
   5 #
   6 # e.g.
   7 #   spring constant k=df/dx      \__
   8 #   velocity        v=dx/dt      /  `--> df/dt = kv --> f = kvt  (+ f0)
   9 #   unfolding rate  K=1 frac/s
  10 #
  11 #   Population as a function of force/time follows
  12 #     p(t) = exp(-tK) = exp(-FK/kv) = p(F)
  13 #   so histogram of unfolding vs. force p(F) normalized to p(0)=1 should follow
  14 #     -binwidth N0 dp/dF = binwidth N0 pK/kv = binwidth N0 K/kv exp(-FK/kv)
  15 #
  16 #   test with lots of runs of
  17 #     sawsim -v1 -mhooke -a1 -kconst -K1 -F1
  18 #   which is not a problem, because 200 runs of that takes ~1 second:
  19 #     time xtimes 200 sawsim -v1 -mhooke -a1 -kconst -K1 -F1 > /dev/null
  20 #
  21 # usage: const_rate.sh <velocity> <spring_constant> <unfolding_rate>
  22
  23 if [ $# -ne 3 ]
  24 then
  25     echo "usage: const_rate.sh <velocity> <spring_constant> <unfolding_rate>"
  26     exit 1
  27 fi
  28
  29 V=$1
  30 K=$2
  31 RATE=$3
  32 DEBUG=0
  33
  34 # since we're on the cluster, we'll go crazy and get really good statistics ;)
  35 N=10000
  36 df=`python -c "print 0.1/$RATE"`
  37
  38 SAWSIM=../../bin/sawsim
  39 DATA=const_rate_${V}_${K}_$RATE.d
  40 HIST=const_rate_${V}_${K}_$RATE.hist
  41 COMPFILE=const_rate_${V}_${K}_$RATE.compare
  42 GPFILE=const_rate_${V}_${K}_$RATE.gp
  43 PNGFILE=const_rate_${V}_${K}_$RATE.png
  44
  45 # set up the theory
  46 # our predicted histogram (see notes above) is
  47 #  hist(F) = binwidth N0 K/kv exp(-FK/kv) = A exp(F/B)
  48 theoryA=`python -c "print $df*$N*$RATE/float($K*$V)"`
  49 theoryB=`python -c "print float(-$K*$V)/$RATE"`
  50 THEORY_FN="$theoryA * exp(x / ($theoryB))"
  51
  52
  53 # run the experiments
  54 > $DATA
  55 if [ "$ISACLUSTER" -eq 1 ]
  56 then
  57     # Sawsim <= 0.5
  58     #qcmd "xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA"
  59     # Sawsim >= 0.6
  60     qcmd xtimes $N $SAWSIM -v$V -s folded,hooke,$K -N1 -s unfolded,null \
  61         -k folded,unfolded,const,$RATE -q folded \
  62         | grep -v '^#' | cut -f1 >> $DATA
  63 else
  64     # Sawsim <= 0.5
  65     #xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA
  66     # Sawsim >= 0.6
  67     xtimes $N $SAWSIM -v$V -s folded,hooke,$K -N1 -s unfolded,null \
  68         -k folded,unfolded,const,$RATE -q folded \
  69         | grep -v '^#' | cut -f1 >> $DATA
  70 fi
  71
  72 # histogram the data
  73 if [ $DEBUG -eq 1 ]; then
  74     echo "cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST"
  75 fi
  76 cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST
  77
  78 # fit the data to an exponential decay
  79 if [ $DEBUG -eq 1 ]; then
  80     echo "cat $HIST | ./fit_exponential"
  81 fi
  82 FIT=`cat $HIST | ./fit_exponential` || exit 1
  83 B=`echo "$FIT" | sed -n 's/time constant:[[:space:]]*//p'`
  84 A=`echo "$FIT" | sed -n 's/initial value:[[:space:]]*//p'`
  85 echo -e "force scale:\t$B\t(expected $theoryB)" > $COMPFILE
  86 echo -e "initial value:\t$A\t(expected $theoryA)" >> $COMPFILE
  87 cat $COMPFILE
  88
  89 # generate gnuplot script
  90 if [ $DEBUG -eq 1 ]; then
  91     echo "generate Gnuplot script"
  92 fi
  93 echo "set terminal png" > $GPFILE
  94 echo "set output '$PNGFILE'" >> $GPFILE
  95 echo "set logscale y" >> $GPFILE
  96 echo "set xlabel 'Force'" >> $GPFILE
  97 echo "set ylabel 'Deaths'" >> $GPFILE
  98 echo "theory(x) = $THEORY_FN" >> $GPFILE
  99 echo "fit(x) = $A * exp(x / $B)" >> $GPFILE
 100 echo "plot  theory(x), '$HIST' with points, fit(x)" >> $GPFILE
 101
 102 if [ $DEBUG -eq 1 ]; then
 103     echo "Generate PNG"
 104 fi
 105 gnuplot $GPFILE