diff --git a/fac/facfile.py b/fac/facfile.py
index f7f77a11e..ebe1813d2 100644
--- a/fac/facfile.py
+++ b/fac/facfile.py
@@ -24,6 +24,12 @@ def rule(self, cmd, inputs, outputs):
             print('<', i, file=self._f)
         for o in outputs:
             print('>', o, file=self._f)
+    def default(self, cmd, inputs, outputs):
+        print('\n|', cmd, file=self._f)
+        for i in inputs:
+            print('<', i, file=self._f)
+        for o in outputs:
+            print('>', o, file=self._f)
     def compile(self, cpp):
         obj = cpp[:-3]+'o'
         print('\n|', config.cxx, config.cxxflags, '-c', '-o', obj, cpp, file=self._f)
@@ -31,25 +37,22 @@ def compile(self, cpp):
         print('>', obj, file=self._f)
     def link(self, maincpp, exe, objects=set([])):
         obj = {maincpp[:-3]+'o'}
+        headers = set([])
         with open(maincpp) as f:
             maincontents = f.read()
             for i in includes_re.findall(maincontents):
                 if i == 'Monte-Carlo/monte-carlo':
                     i = 'utilities'
-                print('# include ', i, file=self._f)
-                if os.path.exists('src/'+i+'.cpp'):
-                    print('# should link with', 'src/'+i+'.cpp', file=self._f)
+                if os.path.exists('src/'+i+'.cpp') or i[-4:] == 'Fast':
                     obj |= {'src/'+i+'.o'}
-                else:
-                    print('# should not link with', 'src/'+i+'.cpp', file=self._f)
-            print('# instructions:', instructions_re.findall(maincontents), file=self._f)
+                    headers |= {'src/'+i+'.h'}
             for i in instructions_re.findall(maincontents):
-                print('# instructions:', repr(i), file=self._f)
                 try:
                     exec(i)
                 except:
-                    print('# instructions failed!', repr(i), sys.exc_info(), file=self._f)
+                    print('# instructions failed!', repr(i), file=self._f)
+                    print('# instructions failed!', sys.exc_info()[1], file=self._f)
         print('\n|', config.cxx, config.linkflags, '-o', exe, ' '.join(obj | objects), file=self._f)
         print('>', exe, file=self._f)
-        for o in obj | objects:
+        for o in obj | objects | headers:
             print('<', o, file=self._f)
diff --git a/fac/figs.py b/fac/figs.py
new file mode 100644
index 000000000..6a6c6e8af
--- /dev/null
+++ b/fac/figs.py
@@ -0,0 +1,92 @@
+#!/usr/bin/python2
+
+import glob, re, string, os, sys
+
+pyfs = sys.argv[1:]
+
+import_re = re.compile(r"^import\s+(\w+)", re.M)
+
+fixed_open = re.compile(r"=\s*open\(['\"]([^'\"]*)['\"]\s*,\s*['\"]w['\"]\s*\)")
+open_changing_output = re.compile(r"open\(['\"]([^'\"]*)['\"]\s*%\s*(\(.*\))(\s*,[\w\s=]+)*\s*,\s*['\"]w['\"]\s*\)")
+fixed_output = re.compile(r"savefig\(['\"]([^'\"]*)['\"](\s*,[\w\s=]+)*\)")
+changing_output = re.compile(r"savefig\(['\"]([^'\"]*)['\"]\s*%\s*(\(.*\))(\s*,[\w\s=]+)*\s*\)")
+arguments = re.compile(r"^#arg\s+(\w+)\s*=\s*(.*)$", re.M)
+
+fixed_open_input = re.compile(r"=\s*open\(['\"]([^'\"]*)['\"]\s*,\s*['\"]r['\"]\s*\)")
+fixed_input = re.compile(r"^[^\n#]*loadtxt\(['\"]([^'\"]*)['\"]\)", re.M)
+changing_loadtxt_noparens = re.compile(r"^[^\n#]*loadtxt\(['\"]([^'\"]*)['\"]\s*%\s*([^\(\)\n]*)(\s*,[\w\s=]+)*\s*\)", re.M)
+changing_loadtxt = re.compile(r"^[^\n#]*loadtxt\(['\"]([^'\"]*)['\"]\s*%\s*(\([^\)]*\))(\s*,[\w\s=]+)*\s*\)", re.M)
+changing_input = re.compile(r"input:\s*['\"]([^'\"]*)['\"]\s*%\s*(\(.*\))")
+list_comprehension_input = re.compile(r"input:\s*(\[.+for.*in.*\])\n")
+
+def friendly_eval(code, context, local = None):
+    try:
+        return eval(code, local)
+    except:
+        print ("\nError evaluating '%s' from file %s" % (code, context))
+        raise
+
+
+for fname in pyfs:
+    f = open(fname, 'r')
+    contents = f.read()
+    f.close()
+
+    inputs = fixed_input.findall(contents)
+    inputs += fixed_open_input.findall(contents)
+
+    imports = import_re.findall(contents)
+    for i in imports:
+        ipath = os.path.join(os.path.dirname(fname), i+'.py')
+        if os.path.exists(ipath):
+            inputs.append(ipath)
+
+    outputs = fixed_output.findall(contents)
+    outputs = [x[0] for x in outputs]
+    outputs += fixed_open.findall(contents) # add any files created with open(...)
+    argvals = arguments.findall(contents)
+    if len(argvals) > 0:
+        coutputs = changing_output.findall(contents) + open_changing_output.findall(contents)
+        list_inputs = list_comprehension_input.findall(contents)
+        cinputs = changing_input.findall(contents)
+        cinputs += changing_loadtxt.findall(contents)
+        cinputs += changing_loadtxt_noparens.findall(contents)
+        allvalues = [{}]
+        commandlineformat = ""
+        commandlineargs = []
+        for arg in argvals:
+            commandlineformat += " \"%s\""
+            commandlineargs.append(arg[0])
+            newallvalues = []
+            for v in friendly_eval(arg[1], fname):
+                for av in allvalues:
+                    newav = av.copy()
+                    newav[arg[0]] = v
+                    newallvalues.append(newav)
+            allvalues = newallvalues
+        for a in allvalues:
+            aa = commandlineformat % friendly_eval('(' + string.join(commandlineargs, ",") + ')', 'file '+ fname, a)
+            extrainputs = []
+            extraoutputs = []
+            for i in list_inputs:
+                #print 'examining input:', i
+                fnames = friendly_eval(i, fname, a)
+                #print 'found', fnames
+                extrainputs += fnames
+            for i in cinputs:
+                extrainputs.append(i[0] % friendly_eval(i[1], fname, a))
+            for o in coutputs:
+                extraoutputs.append(o[0] % friendly_eval(o[1], fname, a))
+            print "? python2 -B %s %s" % (fname, aa)
+            for i in inputs + extrainputs:
+                print '<', i
+            for o in outputs + extraoutputs:
+                print '>', o
+            print
+    else:
+        print '| python2 -B', fname
+        for i in inputs:
+            print '<', i
+        for o in outputs:
+            print '>', o
+        print
diff --git a/fac/latex.py b/fac/latex.py
new file mode 100644
index 000000000..d9bdf9332
--- /dev/null
+++ b/fac/latex.py
@@ -0,0 +1,42 @@
+#!/usr/bin/python2
+
+import re, string, sys, os
+
+import facfile
+
+if len(sys.argv) > 1:
+    texfs = sys.argv[1:]
+else:
+    print("Need arguments")
+    os.exit(1)
+
+documentclassre = re.compile(r'\\documentclass(\[[^\]]+\])?{')
+graphicre = re.compile(r'^\s*\\includegraphics(\[[^\]]+\])?{([^}]+)}', re.MULTILINE)
+inputre = re.compile(r'\\input(\[[^\]]+\])?{([^}]+)}')
+
+for t in texfs:
+    fac = facfile.facfile(os.path.dirname(t)+'/.tex.fac')
+
+    fname = os.path.basename(t)
+    f = open(t, 'r')
+    latex = f.read()
+    inputs = {fname}
+    outputs = {fname[:-3]+'pdf', fname[:-3]+'bbl'}
+
+    os.chdir(os.path.dirname(t)) # so we can look up figures easily
+    def graphics_name(x):
+        if os.path.exists(x):
+            return x
+        if os.path.exists(x+'.png'):
+            return x+'.png'
+        if os.path.exists(x+'.jpg'):
+            return x+'.jpg'
+        if len(x) > 4 and x[-4:] in ['.pdf', '.png', '.jpg']:
+            return x
+        return x+'.pdf'
+    inputs |= set([graphics_name(x[1]) for x in graphicre.findall(latex)])
+    inputs |= set([x[1]+'.tex' for x in inputre.findall(latex)])
+    fac.default('pdflatex %s && pdflatex %s && bibtex %s && pdflatex %s'
+                % (fname, fname, fname[:-4], fname), inputs, outputs)
+
+
diff --git a/fac/monte-carlos.py b/fac/monte-carlos.py
index b42c65cb0..f530339b8 100644
--- a/fac/monte-carlos.py
+++ b/fac/monte-carlos.py
@@ -24,10 +24,10 @@
 
 
 for method in ["nw","bubble_suppression","robustly_optimistic","gaussian",
-               "wang_landau","walker_optimization", 'tmmc']:
+               "wang_landau","walker_optimization", 'tmmc', 'kT 1', 'kT 2']:
     for ww in [1.3, 1.5]:
         for ff in [0.3, 0.8]:
-            for N in range(5,20):
+            for N in range(5,21):
                 outputs = ["papers/square-well-liquid/data/periodic-ww%04.2f-ff%04.2f-N%i-%s-%s.dat"
                            % (ww, ff, N, method.replace(' ',''), postfix)
                            for postfix in ['g', 'E', 'lnw', 's', 'transitions']]
diff --git a/fac/start-here.fac b/fac/start-here.fac
index 9d7402e75..69e00ea9d 100644
--- a/fac/start-here.fac
+++ b/fac/start-here.fac
@@ -17,3 +17,12 @@ C __pycache__
 < config.py
 > ../src/haskell/.haskell.fac
 C __pycache__
+
+| python3 -B latex.py papers/square-well-liquid/histogram-paper.tex
+< config.py
+> ../papers/square-well-liquid/.tex.fac
+C __pycache__
+
+| cd ../papers/square-well-liquid && python2 -B ../../fac/figs.py figs/plot-dos.py figs/plot-histograms.py figs/plot-samples.py figs/plot-scaling.py figs/plot-transitions.py figs/plot-uhc.py > .figs.fac
+> ../papers/square-well-liquid/.figs.fac
+
diff --git a/papers/bilge-latex.py b/papers/bilge-latex.py
deleted file mode 100644
index 46ff0835e..000000000
--- a/papers/bilge-latex.py
+++ /dev/null
@@ -1,26 +0,0 @@
-#!/usr/bin/python2
-
-import glob, re, string, sys
-
-if len(sys.argv) > 1:
-    texfs = sys.argv[1:]
-else:
-    texfs = glob.glob('*.tex')
-
-documentclassre = re.compile(r'\\documentclass(\[[^\]]+\])?{')
-graphicre = re.compile(r'^\s*\\includegraphics(\[[^\]]+\])?{([^}]+)}', re.MULTILINE)
-inputre = re.compile(r'\\input(\[[^\]]+\])?{([^}]+)}')
-
-for fname in texfs:
-    f = open(fname, 'r')
-    latex = f.read()
-    print '| pdflatex %s && pdflatex %s && bibtex %s && pdflatex %s' % (fname, fname, fname[:-4], fname)
-    print '>', fname[:-3]+'bbl'
-    print '>', fname[:-3]+'pdf'
-    if documentclassre.search(latex):
-        for pdf in [x[1]+'.pdf' for x in graphicre.findall(latex)]:
-            print '<', pdf
-        for tex in [x[1]+'.tex' for x in inputre.findall(latex)]:
-            print '<', tex
-    print
-
diff --git a/papers/fuzzy-fmt/figs/radial-wca.cpp b/papers/fuzzy-fmt/figs/radial-wca.cpp
index 9a927ccb9..7391b0164 100644
--- a/papers/fuzzy-fmt/figs/radial-wca.cpp
+++ b/papers/fuzzy-fmt/figs/radial-wca.cpp
@@ -23,7 +23,7 @@
 #include "utilities.h"
 #include "handymath.h"
 
-/*--
+/* --
 
 for kT in np.arange(0.1, 1.0, 0.1):
   for rho in np.arange(0.1, 1.0, 0.1):
@@ -31,7 +31,7 @@ for kT in np.arange(0.1, 1.0, 0.1):
               [exe],
               ["papers/fuzzy-fmt/figs/soft-sphere%06.4f-%04.2f.dat" % (kT, rho)])
 
---*/
+-- */
 
 //const double nm = 18.8972613;
 // Here we set up the lattice.
diff --git a/papers/fuzzy-fmt/figs/soft-sphere.cpp b/papers/fuzzy-fmt/figs/soft-sphere.cpp
index 48ad77dc1..c444bd2a0 100644
--- a/papers/fuzzy-fmt/figs/soft-sphere.cpp
+++ b/papers/fuzzy-fmt/figs/soft-sphere.cpp
@@ -23,7 +23,7 @@
 #include "utilities.h"
 #include "handymath.h"
 
-/*--
+/* --
 
 for kT in np.arange(0.1, 1.0, 0.1):
   for rho in np.arange(0.1, 1.0, 0.1):
@@ -31,7 +31,7 @@ for kT in np.arange(0.1, 1.0, 0.1):
               [exe],
               ["papers/fuzzy-fmt/figs/radial-wca-%06.4f-%04.2f.dat" % (kT, rho)])
 
---*/
+-- */
 
 //const double nm = 18.8972613;
 // Here we set up the lattice.
diff --git a/papers/square-well-liquid/figs/plot-dos.py b/papers/square-well-liquid/figs/plot-dos.py
index 1fbe2b966..5f3a19ae5 100755
--- a/papers/square-well-liquid/figs/plot-dos.py
+++ b/papers/square-well-liquid/figs/plot-dos.py
@@ -20,7 +20,7 @@
 #arg N = [10, 20, 100, 200, 1000]
 
 versions = eval(sys.argv[4])
-#arg versions = [["nw","wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization", "kT2", "kT1"]]
+#arg versions = [["nw","wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization", "kT2", "kT1", "tmmc"]]
 
 # input: ["data/periodic-ww%04.2f-ff%04.2f-N%i-%s-%s.dat" % (ww, ff, N, version, data) for version in versions for data in ["E","lnw"]]
 
diff --git a/papers/square-well-liquid/figs/plot-histograms.py b/papers/square-well-liquid/figs/plot-histograms.py
index e7bb07149..af8b14f9b 100755
--- a/papers/square-well-liquid/figs/plot-histograms.py
+++ b/papers/square-well-liquid/figs/plot-histograms.py
@@ -20,7 +20,7 @@
 #arg N = [10, 20, 100, 200, 1000]
 
 versions = eval(sys.argv[4])
-#arg versions = [["nw","wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization","kT2","kT1"]]
+#arg versions = [["nw","wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization","kT2","kT1","tmmc"]]
 
 # input: ["data/periodic-ww%04.2f-ff%04.2f-N%i-%s-E.dat" % (ww, ff, N, version) for version in versions]
 
diff --git a/papers/square-well-liquid/figs/plot-scaling.py b/papers/square-well-liquid/figs/plot-scaling.py
index 949b78865..e0ff60459 100755
--- a/papers/square-well-liquid/figs/plot-scaling.py
+++ b/papers/square-well-liquid/figs/plot-scaling.py
@@ -6,7 +6,7 @@
 import numpy, glob, re, string
 import styles
 
-if len(sys.argv) not in [4,5]:
+if len(sys.argv) not in [5,6]:
     print 'useage: %s ww ff N versions show' % sys.argv[0]
     exit(1)
 
@@ -16,51 +16,43 @@
 ff = float(sys.argv[2])
 #arg ff = [0.1, 0.2, 0.3, 0.4]
 
-versions = eval(sys.argv[3])
-#arg versions = [["wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization"]]
+all_Ns = eval(sys.argv[3])
+#arg all_Ns = [[5,6,7,8,9,10,20]]
 
-# input: ["data/periodic-ww%04.2f-ff%04.2f-N%i-%s-%s.dat" % (ww, ff, N, version, dat) for version in versions for N in [5,6,7,8,9,10,20] for dat in ['s', 'lnw', 'E']]
+versions = eval(sys.argv[4])
+#arg versions = [["wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization","tmmc"]]
+
+# input: ["data/periodic-ww%04.2f-ff%04.2f-N%i-%s-%s.dat" % (ww, ff, N, version, dat) for version in versions for N in all_Ns for dat in ['s', 'lnw', 'E']]
 
 N_regex = re.compile(r'-N([0-9]+)')
 initialization_iters_regex = re.compile(r'# iterations:\s+([0-9]+)')
 
 plt.title('Scaling for $\lambda=%g$, $\eta=%g$' % (ww, ff))
-Ns = {}
 init_iters = {}
 Emins = {}
 samples = {}
 for version in versions:
-    Ns[version] = []
-    init_iters[version] = []
-    Emins[version] = []
-    samples[version] = []
-    for filename in glob.glob("data/periodic-ww%04.2f-ff%04.2f-N*-%s-lnw.dat" % (ww, ff, version)):
-        N = int(N_regex.findall(filename)[0])
-        Ns[version].append(N)
-        wildfilename = "data/periodic-ww%04.2f-ff%04.2f-N%d-%s-%%s.dat" % (ww, ff, N, version)
-
-        with open(filename, 'r') as content_file:
-            content = content_file.read()
-        init_iters[version].append(int(initialization_iters_regex.findall(content)[0]))
-
-        E_data = numpy.loadtxt(wildfilename % 'E')
-        Emins[version].append(E_data[:, 0].max())
-
-        sample_data = numpy.loadtxt(wildfilename % 's')
-        samples[version].append(sample_data[len(sample_data[:,1])-1, 1])
-
-    # The following sorts the lists according to the number of spheres
-    indexes = range(len(Ns[version]))
-    indexes.sort(key=Ns[version].__getitem__)
-    Ns[version] = map(Ns[version].__getitem__, indexes)
-    init_iters[version] = map(init_iters[version].__getitem__, indexes)
-    Emins[version] = map(Emins[version].__getitem__, indexes)
-    samples[version] = map(samples[version].__getitem__, indexes)
-
-    plt.figure(1)
-    plt.semilogy(Ns[version], init_iters[version], styles.color[version]+'.-', label=version)
-    plt.figure(2)
-    plt.plot(Ns[version], Emins[version], styles.color[version]+'.-', label=version)
+  init_iters[version] = []
+  Emins[version] = []
+  samples[version] = []
+  for N in all_Ns:
+    filename = "data/periodic-ww%04.2f-ff%04.2f-N%d-%s-lnw.dat" % (ww, ff, N, version)
+    wildfilename = "data/periodic-ww%04.2f-ff%04.2f-N%d-%s-%%s.dat" % (ww, ff, N, version)
+
+    with open(filename, 'r') as content_file:
+        content = content_file.read()
+    init_iters[version].append(int(initialization_iters_regex.findall(content)[0]))
+
+    E_data = numpy.loadtxt(wildfilename % 'E')
+    Emins[version].append(E_data[:, 0].max())
+
+    sample_data = numpy.loadtxt(wildfilename % 's')
+    samples[version].append(sample_data[len(sample_data[:,1])-1, 1])
+
+  plt.figure(1)
+  plt.semilogy(all_Ns, init_iters[version], styles.color[version]+'.-', label=version)
+  plt.figure(2)
+  plt.plot(all_Ns, Emins[version], styles.color[version]+'.-', label=version)
 
 plt.figure(1)
 plt.xlabel('$N$')
@@ -88,8 +80,8 @@
 \hline\hline
 """)
 
-for i in range(len(Ns['bubble_suppression'])):
-    N = Ns['bubble_suppression'][i]
+for i in range(len(all_Ns)):
+    N = all_Ns[i]
     tex.write(r""" N = %d \\
   initialization""" % N)
 
diff --git a/papers/square-well-liquid/figs/plot-uhc.py b/papers/square-well-liquid/figs/plot-uhc.py
index 9498f5d49..b2126c655 100755
--- a/papers/square-well-liquid/figs/plot-uhc.py
+++ b/papers/square-well-liquid/figs/plot-uhc.py
@@ -22,7 +22,7 @@
 #arg N = [20, 100, 200, 1000]
 
 versions = eval(sys.argv[4])
-#arg versions = [["nw","wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization","kT2","kT1"]]
+#arg versions = [["nw","wang_landau","robustly_optimistic","gaussian","bubble_suppression","walker_optimization","kT2","kT1","tmmc"]]
 
 # input: ["data/periodic-ww%04.2f-ff%04.2f-N%i-%s-%s.dat" % (ww, ff, N, version, data) for version in versions for data in ["E","lnw"]]
 
diff --git a/papers/square-well-liquid/figs/styles.py b/papers/square-well-liquid/figs/styles.py
index e31284782..5695ff828 100644
--- a/papers/square-well-liquid/figs/styles.py
+++ b/papers/square-well-liquid/figs/styles.py
@@ -2,6 +2,7 @@
           'kT1': 'b',
           'kT2': 'c',
           'kT0.1': 'g',
+          'tmmc': 'k',
           'bubble_suppression': 'k',
           'wang_landau': 'g',
           'vanilla_wang_landau': 'b',
@@ -13,6 +14,7 @@
          'kT1': '--',
          'kT2': ':',
          'kT0.1': '-.',
+         'tmmc': '-',
          'bubble_suppression': '-',
          'wang_landau': '-',
          'vanilla_wang_landau': '-',
@@ -24,6 +26,7 @@
           'kT1': '$kT/\epsilon = 1$ sim.',
           'kT2': '$kT/\epsilon = 2$ sim.',
           'kT0.1': '$kT/\epsilon = 0.1$ sim.',
+          'tmmc': 'tmmc',
           'bubble_suppression': 'bubble suppression',
           'wang_landau': 'Wang-Landau',
           'vanilla_wang_landau': 'Vanilla Wang-Landau',
diff --git a/src/Monte-Carlo/square-well.cpp b/src/Monte-Carlo/square-well.cpp
index 138fd15cb..b5992234e 100644
--- a/src/Monte-Carlo/square-well.cpp
+++ b/src/Monte-Carlo/square-well.cpp
@@ -782,11 +782,17 @@ void sw_simulation::update_weights_using_transitions(double fractional_precision
     // compute D_n = T*D_{n-1}
     dos_magnitude_squared = 0;
     for (int i = 0; i < energies_observed; i++){
-      double TD_i = 0;
+      double TD_i = 0, norm_i = 0, dos_min = 1e300;
       for(int e = max(0,i-biggest_energy_transition);
-          e <= min(energy_levels,i+biggest_energy_transition); e++)
-        TD_i += transitions(i,e-i)*exp(ln_old_dos[e]);
-      ln_dos[i] = log(TD_i);
+          e <= min(energy_levels,i+biggest_energy_transition); e++) {
+        if (ln_old_dos[e] < dos_min) dos_min = ln_old_dos[e];
+      }
+      for(int e = max(0,i-biggest_energy_transition);
+          e <= min(energy_levels,i+biggest_energy_transition); e++) {
+        TD_i += transitions(i,e-i)*exp(ln_old_dos[e] - dos_min);
+        norm_i += transitions(i,e-i);
+      }
+      ln_dos[i] = log(TD_i/norm_i) + dos_min;
       dos_magnitude_squared += exp(2*ln_dos[i]);
     }