Saturday, 16 February 2013

Simple Glossaries in LaTeX

I use LaTeX (via MikTeX) and TeXnicCenter a bit, starting from when I did my PhD. It's a great document preparation system which has some annoying sides, but I won't go down that path. I tried to set up a glossary in a document recently, and it wasn't obvious to me how to do that. I did get it working though. Here's what I did, so that others (including myself) don't have to figure it out again:
  1. Follow the steps provided at LaTeX Matters, namely:
    1.  Put two lines into your preamble (before \begin{document}):
    2. When you define or introduce certain terms, wrap them in a \glossary command like so:
      \glossary{name={entry name}, description={entry description}}
      The inside curly braces are only necessary if you have commas in the entries.
    3. Place this command where you want the glossary to appear:
    4. If you want your glossary to have an entry in the Table of Contents, then put this with the previous line:
    5. Then, instead of running a command as described in the post, follow the next step, which is equivalent, but for TeXnicCenter.
  2. Prompted by Archiddeon's comment in a LaTex forum, in TeXnicCenter, add a post-processing step to a build profile:
    1. Select the Build menu then choose Define Output Profiles...
    2. On the left, select a profile you want to start from
    3. Click the Copy button down below, enter a new name for the profile and click OK
    4. Then click on the Postprocessor tab
    5. Create a new postprocessor step by clicking the little "new folder" icon at the top right
    6. Then enter a name like "Prepare the glossary" and the following:
      Executable: makeindex
      Arguments: "%tm.glo" -s "" -t "%tm.glg" -o "%tm.gls"
  3. Click OK until you're back at the main window!
Then you should be able to select the new profile and generate the document with a glossary! I'd run it 2-3 times to let the page numbers settle. Joy!

Friday, 8 February 2013

The 147 Problem from Programming Praxis

Programming Praxis is a great site. The 147 problem was recently posted there, and I tackled it with a bit too much time... Following is the code, with discussion.

(defun sum1/ (s)
    "Read the name as Sum the Inverses"
    (apply #'+ (mapcar #'/ s)))

(defun limits (n)
    "I took the limits discussion on the site a bit further..."
    (loop for i from n downto 1
          for prod = (apply #'* prods); Produces 1 to start with
          for prods = (cons (1+ prod) prods)
          for min = (min (car prods) n)
          for max = (max (* i prod) n)
        collect (list min max)))

When calling limits with n = 5, we get ((2 5) (3 8) (5 18) (5 84) (5 1806)) which shows the maximum ranges of each number we'll search for. The actual search for solutions to 1/a+1/b+1/c+1/d+1/e = 1 then looks like this:

(defun search1 (n &optional givens limits)
    "Actually finds all of the valid solutions to the 147
problem with n numbers. Call it without the optional arguments
to start it off."
    (let ((limits (if (null givens) (limits n) limits)))

        ;; Limits is whittled down during the recursion,
        ;; so detect starting off by looking at givens.
        (if (= n (length givens))
            (if (= 1 (sum1/ givens))
                (list (reverse givens)))
            (loop for next from (max (caar limits)

                                     (if givens (car givens) 0))
                           to (cadar limits)
                  for s = (cons next givens)
                until (or (every (lambda (i) (> i n)) s)
                          (and (not (cdr limits)) (< (sum1/ s) 1)))
                appending (search1 n s (cdr limits))))))

This code times how long it takes to solve the first few solutions before the times get too big:

(time (loop for n from 0 to 5 collect (nreverse (search1 n))))

Which gave 4.1 seconds on my computer.

Note that the code works for n = 0 and 1 (producing the empty list, and a single solution: 1, respectively). As some other people have mentioned, Common Lisp is well thought out. I think this also supports that thought, as I didn't design the algorithm specifically to handle cases n = 0 and 1.

And lastly, the following code finds the smallest distinct solutions, as described at Programming Praxis:

  (remove nil
    (mapcar (lambda (s)
              (if (apply #'/= s); Distinctness
                (cons s (apply #'+ s)))); Sum the denominators
            (search1 5)))
  #'< :key #'cdr); Sort by the sum of denominators

Which gives (3 4 5 6 20), confirming our solution.

Happy coding!

Saturday, 2 February 2013

Solving the Discretised Traffic Model Equations

This post follows on from when we discretised our traffic model equations. We now solve them in this post.

The solving process is quite simple since we avoided most references to the future values of variables. If we did, we'd have to use matrices to solve the discretised model equations. As it is, we can simply calculate each value separately.

The code is modeled around each cell, so here's the discretised cell equation: $$\rho(x_i,t_f) = \rho(x_i,t_c) - \frac{t_f-t_c}{x_{i+}-x_{i-}}\left(\rho(x_{i+},t_c) v(x_{i+},t_c) - \rho(x_{i-},t_c) v(x_{i-},t_c)\right).$$ where
\(\rho\) = vehicular density (veh/m);
\(t\) = time (s);
\(t_c\) = the current time (s);
\(t_f\) = the future time (s);
\(x\) = distance (m) along road;
\(i\) = the index of the current cell, ranging from 1 to \(n\);
\(x_i\) = the location of the 'centroid' of the current cell (m);
\(x_{i-}\) = the left (upstream) boundary of the current cell
\(x_{i+}\) =
\(v_f\) = free flow speed (m/s);
\(\rho_j\) = jam density (veh/m), the average vehicles per metre in stationary traffic; and
\(v\) = vehicle speed (m/s), and is given by: \(v = v_f(1-\rho/\rho_j)\).

The left boundary is given by: $$\rho(0,t) = \rho_j/4.$$

Oh, and in case it isn't clear, the flow of traffic is in the direction of increasing \(x\) , which is to the right. That means that upstream is to the left, and downstream is to the right.

The above equations have been implemented in Common Lisp, with MathP (but I'll wait till the end to show the final code). It outputs the values scaled up by 1000 for convenience. To call it, we evaluate something like this:

(step-simulation :road-length 1000 :n 20 :dt 1.0 :time-run 200)

Great! So, let's look at the results. Here is a visualisation of the results for the above simulation:

Initial traffic simulation results. Time increases from top to bottom, distance increases from left to right. Light traffic is green, medium traffic is yellow, and heavier traffic is red.

Not so good... there are fluctuations in the solution, causing negative vehicle densities. Thankfully, that can be fixed with an Upwind scheme. In our case, we just move the centroid more upstream within the cell. This removes the fluctuations, but probably increases error elsewhere in the simulation. The results with the Upwind scheme are:
"Upwind" traffic simulation results. Time increases from top to bottom, distance increases from left to right. Light traffic is green, medium traffic is yellow, and heavier traffic is red.

Looking good. The traffic on the left fills the road over time, and there are no fluctuations.

The final code to implement this is:

(defparameter *free-flow-speed* (/ 60 3.6); = 60 km/h
    "Speed (m/s) of freely flowing, low density traffic.")
(defparameter *jam-density* (/ 7.0); = 1 veh every 7 m

    "Jam density (veh/m), the average vehicles per metre in stationary traffic.")
(defparameter *road-length* 1000.0; = 1 km

    "Length (m) of road being simulated.")
(defparameter *cell-count* 10

    "Number of cells road is split into for simulation.")

(defun step-simulation (&key
        (vf *free-flow-speed*)
        (pj *jam-density*)
        (road-length *road-length*)
        (n *cell-count*)
        (dt 1.0); Time step in seconds.
        (time-run 10.0); Time to run simulation for.
        (t0 0.0); Starting time for simulation.
        (upwind 0.5); Upwinding parameter [0 = only upstream, 1 = only downstream]
        p0; Initial road density for each cell.
        pin); Input density function (veh/m). Should take one argument - the time.
    dx = roadLength/(n-1)
; Spacing between cell centres.
    pin = if pin pin fn(_) pj/4.0

    ;; Set the boundary density in the initial state.
    unless(p0 p0 := makeArray(n :initialElement 0.0) p0[0] := #!pin(t0))
    v(p) = vf*(1-p/pj); Velocity as function of density.
    x(i) = (i-1)*dx; Cell centroid locations.
    xp(i) = if i<n (i-0.5)*dx roadLength; Plus, or right (downstream) boundary locations.
    xm(i) = if i>1 (i-1.5)*dx 0; Minus, or left (upstream) boundary locations.
    stepsim(pold atime dt) =
        makeVector(n fn(i) {
            pm = if i>1 (1-upwind)*pold[i-2]+upwind*pold[i-1] pold[0]
            pp = if i<n (1-upwind)*pold[i-1]+upwind*pold[i] pold[n-1]
            in = max(0 dt*pm*v(pm))
            out= max(0 dt*pp*v(pp))
            if i==1 #!pin(atime) pold[i-1] + (in - out)/(xp(i)-xm(i))})
(loop for atime from t0 to (+ t0 time-run) by dt
            for ps = p0 then (stepsim ps atime dt)
        do (loop for p across ps
            do (format t " ~4,1,3F" p)
            finally (terpri)))})