1

# dxo: A System for Relational Algebra and Differentiation

Julie S. Steele  Georgetown Day School, USA<EMAIL_ADDRESS>and
William E. Byrd  University of Alabama at Birmingham, USA<EMAIL_ADDRESS>
###### Abstract.

We present _dxo_ , a relational system for algebra and differentiation,
written in miniKanren. _dxo_ operates over math expressions, represented as
s-expressions. _dxo_ supports addition, multiplication, exponentiation,
variables (represented as tagged symbols), and natural numbers (represented as
little-endian binary lists). We show the full code for _dxo_ , and describe in
detail the four main relations that compose _dxo_. We present example problems
_dxo_ can solve by combining the main relations. Our differentiation relation,
do, can differentiate polynomials, and by running backwards, can also
integrate. Similarly, our simplification relation, simpo, can simplify
expressions that include addition, multiplication, exponentiation, variables,
and natural numbers, and by running backwards, can complicate any expression
in simplified form. Our evaluation relation, evalo, takes the same types of
expressions as simpo, along with an environment associating variables with
natural numbers. By evaluating the expression with respect to the environment,
evalo can produce a natural number; by running backwards, evalo can generate
expressions (or the associated environments) that evaluate to a given value.
reordero also takes the same types of expressions as simpo, and relates
reordered expressions.

relational programming, differentiation, simplification, miniKanren, Racket,
Scheme

††copyright: none††price: ††ccs: Computing methodologies Computer algebra
systems††ccs: Software and its engineering Functional languages††ccs: Software
and its engineering Constraint and logic languages

##  1\. Introduction

Consider this calculus problem:

> Find two different polynomials, $f(x)$ and $g(x)$, and two different natural
> numbers $a$ and $b$, such that $f^{\prime}(a)=b$, and $g^{\prime}(b)=a$.

Differentiating polynomials is an easy calculus problem, but the problem above
is more complicated because of the relationships between the polynomials,
their derivatives, and the two natural numbers. We invite the reader to pause,
try to find solutions to this problem, and to think about how these types of
problems might be solved more generally.

We have developed a relational algebra system, _dxo_ , that uses relational
programming to solve problems like the one above. We show the `run` expression
for solving this problem in Section 2. _dxo_ is a collection of four main
relations: `simpo` for simplification, `do` for differentiation, `evalo` for
evaluation, and `reordero` for permuting arguments. Implementing _dxo_
relationally makes it flexible. For example, the relation `do` can
differentiate polynomials with respect to some variable. Since `do` is a
relation, it can also integrate polynomials. Also, the expression to be
differentiated and its derivative can both contain fresh logic variables. The
relations `simpo`, `evalo`, and `reordero` similarly benefit from this
flexibility.

We assume the reader is familiar with core miniKanren (Friedman et al., 2018;
Byrd, 2009; Byrd and Friedman, 2006) (`==`, `fresh`, `conde`, `run`), extended
with disequality (`=/=`) and `absento` constraints (Byrd et al., 2012).
Detailed explanations to the core miniKanren language can be found in Friedman
et al. (2018), Byrd (2009), and Byrd and Friedman (2006). Descriptions of
disequality and `absento` constraints can be found in Byrd et al. (2012) and
Byrd et al. (2017).

Section 2 gives a high-level explanation of _dxo_ , its uses, and its four
main relations. Section 3 explains in detail the main relations. Section 4
discusses some open problems and possible future work. Section 5 discusses
related work. We conclude the paper in Section 6. Appendix A contains the full
implementation of _dxo_.

##  2\. High-Level Overview

_dxo_ is composed of four main relations, `simpo`, `do`, `evalo`, and
`reordero`, that when used in combination can solve interesting
differentiation math problems. Here are the four relations and their uses:

(simpo comp simp):

relates `comp` and `simp`, where `comp` can be any arithmetic expression and
`simp` is an equivalent, fully simplified one;

(do x expr deriv):

relates a polynomial expression `expr` with its derivative `deriv`, where the
derivative is with respect to `x`;

(evalo env expr value):

relates an expression `expr` with its value `value`, where each variable in
`expr` is associated with a natural number by the environment `env`;

(reordero e1 e2):

relates two equivalent expressions, `e1` and `e2`, by changing the order of
subexpressions in an addition or multiplication in any level of the other
expression.

Figure 1 contains the grammar for expressions accepted by `simpo`, `evalo`,
and `reordero`, and Figure 2 contains the grammar for polynomial expressions
accepted as the `expr` for `do`. `deriv` is a subset. The implementation of
_dxo_ uses the relational arithmetic system created by Oleg Kiselyov, which is
presented in Friedman et al. (2018) and Kiselyov et al. (2008).

¡dxo-expression¿ ::= ¡numeral-or-variable¿ ‘(+ ’ ¡dxo-expression¿ $\ldots$ ‘)’
‘(* ’ ¡dxo-expression¿ $\ldots$ ‘)’ ‘(^’ ¡dxo-expression¿ ¡dxo-expression¿ ‘)’

¡numeral-or-variable¿ :: = ¡tagged-numeral¿ — ¡tagged-variable¿

¡tagged-variable¿ ::= ‘(var ’ ¡symbol¿ ‘)’

¡tagged-numeral¿ ::= ‘(num ’ ¡numeral¿ ‘)’

¡numeral¿ ::= ‘()’ — ‘(0 . ’ ¡positive-numeral¿ ‘)’ — ‘(1 . ’ ¡numeral¿ ‘)’

¡positive-numeral¿ ::= ‘(0 . ’ ¡positive-numeral¿ ‘)’ — ‘(1 . ’ ¡numeral¿ ‘)’

Figure 1. Grammar for general _dxo_ expressions accepted by simpo, evalo, and
reordero.

¡polynomial-expression¿ ::= ¡numeral-or-variable¿ ‘(+ ’ ¡polynomial-
expression¿ $\ldots$ ‘)’ ‘(* ’ ¡polynomial-expression¿ $\ldots$ ‘)’ ‘(^’
¡numeral-or-variable¿ ¡tagged-numeral¿ ‘)’

Figure 2. Restricted grammar for polynomial expressions accepted by do.

Using the _dxo_ relations, we can solve the problem proposed in the
introduction: find two different polynomials, $f(x)$ and $g(x)$, and two
different natural numbers $a$ and $b$, such that $f^{\prime}(a)=b$, and
$g^{\prime}(b)=a$. We relate `f` and `g` with their derivatives, `fd` and
`gd`, using `do`. Then we use `evalo` to evaluate these derivatives at `a` and
`b` respectively (we do this by making one environment where `x` is `a` and
one where `x` is `b`), and set the evaluation to `b` and `a` respectively.
Last, we make sure `f` and `g` are different but both simplified and `a` and
`b` are different.

    (run 20 (f g envb enva)
      (fresh (b a gd fd)
        (=/= f g)
        (=/= b a)
        (== ‘((x . ,b)) envb)
        (== ‘((x . ,a)) enva)
        (do ’x f fd)
        (simpo f f)
        (do ’x g gd)
        (simpo g g)
        (evalo enva fd b)
        (evalo envb gd a)))
    $\Rightarrow$
    ’(
      ...
      ((num _.0) (var x)                         ; f$=c\textrm{ (where }c\textrm{ is any natural number),}$ g$=x$
       ((x)) ((x 1)))                            ; b$=0$, a$=1$
      ...
      ((var x) (^ (var x) (num (0 0 1 1)))       ; f$=x$, g$=x^{12}$
       ((x 1)) ((x 0 0 1 1)))                    ; b$=1$, a$=12$
      ...
    )

Of the 20 outputs produced by the run expression, many had $b=0$ so we only
showed two. The first shown answer shows

$\frac{d}{dx}[c]=0,\textrm{where }c\textrm{ is any natural number}\textrm{ and
}\frac{d}{dx}[x]=1,\textrm{where }x=5$

The second shown answer shows

$\frac{d}{dx}[x]=1,\textrm{where }x=12\textrm{ and
}\frac{d}{dx}[x^{12}]=12,\textrm{where }x=1$

The concise `run` expression solving this problem shows how _dxo_ benefits
from the expressiveness of relational programming.

We showed a combination of the four main _dxo_ relations in solving the
problem in the introduction. We will shortly demonstrate another way to
combine these core relations in the definition of `anydo`, below.

Let’s use `do` to differentiate the polynomial $x^{3}+x^{0}$. Mathematically,

$\frac{d}{dx}[x^{3}+x^{0}]=(x^{2}*3)+0$

The equivalent call to `do` succeeds:

    (do ’x
        ’(+ (^ (var x) (num (1 1))) (^ (var x) (num ())))        ; $x^{3}+x^{0}=$ expr
        ’(+ (* (^ (var x) (num (0 1))) (num (1 1))) (num ())))   ; $(x^{2}*3)+0=$ deriv

The derivative $(x^{2}*3)+0$ is equivalent to $1*x^{2}*3$, so we might expect
the call

    (do ’x
        ’(+ (^ (var x) (num (1 1))) (^ (var x) (num ())))      ; $x^{3}+x^{0}=$ expr
        ’(* (num (1)) (^ (var x) (num (0 1))) (num (1 1))))    ; $1*x^{2}*3=$ deriv

to succeed. Unfortunately, this call fails because `do` requires the
derivative to be in canonical form, $(x^{2}*3)+0$ in this case. This means
some mathematically correct expression and derivative pairs fail as arguments
to `do`.

We created `anydo` to fix this problem. `(anydo expr deriv x)`, like `do`,
relates an expression with its derivative with respect to `x`, except `anydo`
generalizes this to simplified, complicated, or reordered forms of `expr` and
`deriv`. This relaxes the restriction on `deriv` being in canonical form,
making running “backward” more convenient. Calling `anydo` with the same
arguments as above succeeds:

    (anydo ’(+ (^ (var x) (num (1 1))) (^ (var x) (num ())))     ; $x^{3}+x^{0}=$ expr
           ’(* (num (1)) (^ (var x) (num (0 1))) (num (1 1)))    ; $1*x^{2}*3=$ deriv
           ’x)

`anydo` is centered around a call to `do` with arguments similar to `expr` and
`deriv`, `ecomp` and `dcomp`. `expr` and `ecomp` are similar in that they
simplify to the same value, `esimp`, making them equivalent. `anydo` does the
same for `deriv` and `dcorder`, with the additional step of reordering
`dcorder` to be `dcomp`.

    (define anydo
      (lambda (expr deriv x)
        (fresh (esimp dsimp ecomp dcomp dcorder)
          (simpo expr esimp)
          (simpo ecomp esimp)
          (do x ecomp dcomp)
          (reordero dcomp dcorder)
          (simpo dcorder dsimp)
          (simpo deriv dsimp))))

##  3\. _dxo_ Implementation Walk-through

In this section we explain in detail the four main relations in _dxo_.111We
have released the _dxo_ code under an MIT licence at
https://github.com/JShermanSteele/dxo .

###  3.1. simpo

`(simpo comp simp)` relates `comp` and `simp`, where `comp` can be any
arithmetic expression and `simp` is an equivalent, fully simplified one.
_Simplified_ means making all following simplifications:

  * •

$v+0=v$;

  * •

$v*0=0$;

  * •

$v*1=v$;

  * •

$v^{0}=1$ ($v\neq 0$);

  * •

$v^{1}=v$;

  * •

$0^{v}=0$ ($v\neq 0$);

  * •

and $1^{v}=1$;

where $v$ is any expression. For example, let’s simplify $0^{5}+(2*1)$:

    (run* (simp) (simpo ‘(+ (^ (num ()) (num (1 0 1) ))             ; $0^{5}+2*1=$ comp
                            (* (num (0 1)) (num (1))))
                        simp))
    $\Rightarrow$
    ’((num (0 1)))                                                  ; $2=$ simp

`simpo` has base cases of `(== comp simp)` for `comp` and `simp` being the
same number or variable. The three non-base cases, addition, multiplication,
and exponentiation, deeply recursively simplify sub-expressions by checking
for those simplifyable cases.

If `comp` is ground, `(simpo comp simp)` will either succeed exactly once or
fail because there is at most one way to simplify any concrete expression, and
`simpo` has no overlapping cases when running ”forwards”. If `comp` is fresh,
then `simpo` could succeed, but if it is an impossible relation, simpo can try
longer and longer `comp`s, never succeeding, and loop forever.

An example of these behaviors is that running `(simpo comp simp)` with `comp`
as $1^{1}$ and `simp` as a logic variable succeeds because $1^{1}$ simplified
is $1$. Running with `comp` as a logic variable and `simp` as (the
unsimplified) $1^{1}$ diverges, searching for a `comp` forever.

    (run* (simp) (simpo ’(^ (num (1)) (num (1))) simp))       ; $1^{1}=$ comp
    $\Rightarrow$
    ’((num (1)))

    (run 1 (comp) (simpo comp ’(^ (num (1)) (num (1)))))       ; $1^{1}=$ simp

We can also construct a case with a ground `comp` and partially ground `simp`
in which `simpo` will fail. If `comp` is $1^{1}$, which simplifies to $1$, and
`simp` is any addition expression, which may include fresh variables, `simpo`
will fail finitely.

    (run* (q) (simpo ‘(^ (num  (1)) (num (1)))             ; $1^{1}=$ comp
                     ‘(+ . ,q)))                           ; simp is some addition expression
    $\Rightarrow$
    ’()

###  3.2. do

`(do x expr deriv)` relates a polynomial expression `expr` with its derivative
`deriv`, with respect to `x`. For example, running `do` with `expr` and
`deriv` fresh finds integral/derivative pairs:

    (run 24 (expr deriv) (do ’x expr deriv))

    $\Rightarrow$

    ‘(...
      ((^ (var x) (num (0 1)))                          ; $x^{2}=$ expr
       (* (^ (var x) (num (1))) (num (0 1))))           ; $x^{1}*2=$ deriv
      ...
      ((^ (var x) (num (1 _.0 . _.1)))                    ; $x^{a}\textrm{where }a\textrm{ is odd}=$ expr
       (*                                                 ; $x^{a-1}*a=$ deriv
         (^ (var x) (num (0 _.0 . _.1)))
         (num (1 _.0 . _.1))))
      ...
      ((* (^ (var x) (num ())) (^ (var x) (num ())))      ; $x^{0}*x^{0}=$ expr
       (+                                                 ; $0*x^{0}+x^{0}*0=$ deriv
         (* (num ()) (^ (var x) (num ())))
         (* (^ (var x) (num ())) (num ()))))
    )

The best way to understand `do` is as a case analysis on `expr`, which is
either a variable, a number, an exponentiation, an addition, or a
multiplication.

Since the derivative of a sum is the sum of the derivatives of its parts, when
`expr` and `deriv` are sums, the sub-expressions of `expr` and `deriv` are
pair-wise related using `do`. Since the sum can have any positive number of
terms, a helper relation, `map-do-o`, relates each pair in the sums.

If `expr` is a multiplication, `do` must use the multiplication rule that

$\frac{d}{dx}(ab)=\frac{da}{dx}*b+a*\frac{db}{dx},$

and recur down the list of sub-expressions being multiplied. To improve the
efficiency this process, we wrote a helper relation, `multruleo`, that relates
the list of sub-expressions being multiplied and the multiplication’s
derivative. If the list has length greater than one, `multruleo` separates the
first term $e1$ from the rest, $e2*e3*\ldots$. Applying the multiplication
rule to $e1*e2*e3*\ldots$ yields
$\frac{d}{dx}[e1]*e2*e3*\ldots+e1*\frac{d}{dx}[e2*e3*\ldots]$, which is
recursive with `multruleo` because $\frac{d}{dx}[e2*e3*\ldots]$ is the related
derivative argument to `multruleo` with $e2*e3*\ldots$ as the first argument.
This is `conde` the clause in `do` for multiplication.

    ((fresh (l e)
       (== expr ‘(* . ,l))
       (letrec ((multruleo
                 (lambda (l dd)
                   (fresh (e e* d d* a b)
                     (conde
                       [(== l ‘(,e))(do x e dd)]
                       [(== l ‘(,e . ,e*))
                        (== e* ‘(,a . ,b))
                        (== dd ‘(+ (* ,d . ,e*) (* ,e ,d*)))
                        (do x e d)
                        (multruleo e* d*)])))))
         (multruleo l deriv))))

We could recur through the multiplication by recurring with every shorter
multiplication as an argument to `do`, but our approach is simpler because it
does not exit `multruleo` while reccuring through `expr`’s multiplication.

If `expr` is an exponentiation, the second subexpression in the exponent must
be a number by `do`’s grammar, so $\frac{d}{dx}[x^{n}]=(n*(x^{n-1}))$ where
$n$ is any number. There are three clauses for constants,
$\frac{d}{dx}[x^{0}]=0$, $\frac{d}{dx}[n^{m}]=0$, and $\frac{d}{dx}[n]=0$
where $n$ and $m$ are any numbers so long as they both are not $0$. Finally,
the derivative of just `x` is one.

Since `do` orders `deriv` a certain way (for example, $x^{2}*6$ instead of
$6*x^{2}$), some integratable `deriv`s will fail. This is why `do` should be
used with `reordero`. For example,

    (run 1 (expr)
      (do ’x expr ’(* (^ (var x) (num (1))) (num (0 1))))) ; $x^{1}*2=$ deriv
    $\Rightarrow$
    ’((^ (var x) (num (0 1))))                             ; $x^{2}=$ expr

produces an answer, but switching `deriv`’s multiplication order diverges:

    (run 1 (expr)
      (do ’x expr ’(* (num (0 1)) (^ (var x) (num (1)))))) ; $2*x^{1}=$ deriv

Similarly to `simpo`, `do` succeeds exactly once or fails running “forwards”
when `expr` and `x` are ground. With `expr` fresh, `do` can succeed or can
loop infinitely just like `simpo`.

###  3.3. evalo

`evalo` evaluates, and is useful for solving equations.`(evalo env expr
value)` relates an expression `expr` with its value `value`, where each
variable in `expr` is associated with a natural number by the environment
`env`. For example we can look for expressions that evaluate to $8$ in an
environment that binds `z` to $2$:

    (run 200 (expr) (evalo ‘((z . (0 1))) expr ’(0 0 0 1)))      ; z$=2,$value$=8$
    $\Rightarrow$
    ‘(...
      (* (var z) (num (0 0 1)))                                    ; $z*4=$ expr
      ...
      (^ (num (0 1)) (num (1 1)))                                  ; $2^{3}=$ expr
      ...
      (+ (var z) (num ()) (num (0 1 1)))                           ; $z+0+6=$ expr
    )

`evalo` deeply recurs through `expr`, evaluating tagged little-endian binary
lists into miniKanren numbers. For `evalo`’s base cases when `expr` is a
variable or number, `value` is the variable’s value from `env` or the
miniKanren number respectively.If `expr` is an addition, multiplication, or an
exponentiation, then `evalo` relates the first term with its value `evc`, the
rest with its value `evrest`, and adds, multiplies, or exponentiates `evc` and
`evrest` to obtain `value`.The `evalo` code for addition does this:

    ((== ‘(+ ,c . ,rest) expr)
     (conde
       ((== ’() rest) (evalo env c value))
       ((=/= ’() rest)
       (evalo env c evc)
       (evalo env ‘(+ . ,rest) evrest)
       (pluso evc evrest value))))

This will recur through `rest` and sum all the parts to relate `expr` with
`value`.

An interesting use of `evalo` is to solve algebra problems by making `env`
fresh, for example looking for Pythagorean triples. So we set up
$x^{2}+y^{2}=z^{2}$ and also a $z*z=\mathit{z\mathchar 45\relax squared}$
relation to make sure that `z` is a natural number to find the classic
Pythagorean triple!

    (run 1 (env)
      (fresh (xv yv zv z2v)
        (== ‘((x . ,xv) (y . ,yv) (z . ,zv) (z-squared . ,z2v)) env)
        (evalo env ‘(+ (^ (var x) (num (0 1))) (^ (var y) (num (0 1)))) z2v)
        (*o zv zv z2v)))
    $\Rightarrow$
    ’(((x) (y) (z) (z-squared)))                                           ; x$=0,$y$=0,$z$=0$

Not what we wanted, alas. Setting non-zero constraints, though, produces:

    (run 1 (env)
      (fresh (xv yv zv z2v)
        (poso xv)
        (poso yv)
        (poso zv)
        (== ‘((x . ,xv) (y . ,yv) (z . ,zv) (z-squared . ,z2v)) env)
        (evalo env ‘(+ (^ (var x) (num (0 1))) (^ (var y) (num (0 1)))) z2v)
        (*o zv zv z2v)))
    $\Rightarrow$
    ’(((x 1 1)                                                             ; x$=3$
       (y 0 0 1)                                                           ; y$=4$
       (z 1 0 1)                                                           ; z$=5$
       (z-squared 1 0 0 1 1)))                                             ; z${}^{2}=25$

which is a 3-4-5 right triangle.

If either `env` or `expr` is fresh, `evalo` can loop forever, trying more and
more complicated `env`s or `expr`s. If both `env` or `expr` are ground, then
`evalo` will terminate since `evalo` runs simply forwards in this case.

###  3.4. reordero

`(reordero e1 e2)` relates two equivalent expressions, `e1` and `e2`, by
changing the order of subexpressions in an addition or multiplication in any
level of the other expression. It is useful for taking integrals with `do`. We
can use `reordero` to find all reorderings of an expression:

    (run* (e2) (reordero ‘(+ (num (1)) (* (num (0 1)) (num (1 1)))) e2))   ; $1+2*3=$ e1
    $\Rightarrow$
    ’((+ (num (1)) (* (num (0 1)) (num (1 1))))                            ; $1+2*3=$ e2
      (+ (* (num (0 1)) (num (1 1))) (num (1)))                            ; $2*3+1=$ e2
      (+ (num (1)) (* (num (1 1)) (num (0 1))))                            ; $1+3*2=$ e2
      (+ (* (num (1 1)) (num (0 1))) (num (1))))                           ; $3*2+1=$ e2

`(reordero e1 e2)` relates `e1` and `e2` by having the same outer operation
(addition, multiplication, or exponentiation). For addition and multiplication
the code is,

    ((fresh (o e1* e2*)
       (== ‘(,o . ,e1*) e1)
       (== ‘(,o . ,e2*) e2)
       (typeo o ’+or*)
       (reorderitemso e1* e2*)))

where `e1*` and `e2*` are permutations of each other. `reorderitemso` checks
that `e1*` and `e2*` have the same length, and then calls `reorderinnero` on
them. We created `reorderitemso` to improve speed and divergence behavior of
`reordero` by requiring `e1*` and `e2*` have the same length before
considering the relations in `reorderinnero`. `reorderinnero` relates permuted
lists at any depth. To deeply reorder, `reorderinnero` calls `reordero` on the
corresponding sub-expressions for `e1*` and `e2*`.

    (define reorderinnero
      (lambda (e1* e2*)
        (fresh (c1 rc1 rest1 rest2)
          (conde
            ((== ’() e1*) (== ’() e2*))
            ((== ‘(,c1 . ,rest1) e1*)
             (removeo rc1 e2* rest2)
             (reorderinnero rest1 rest2)
             (reordero c1 rc1))))))

`reordero` greatly reduces infinite loops because `reorderitemso` checks that
its arguments are the same length. This check keeps `e1` and `e2` the same
structure and length at every depth, keeping the search finite. Checks like
this would be useful to add other places in _dxo_ to reduce divergence.

##  4\. Open Problems

`dxo` could be improved by expanding the grammars, improving the speed and
termination, and using automatic differentiation. We would like to add
expressions like $2^{x}$, $sin(x)$, and multiple variables. Currently, _dxo_
searches inefficiently, especially combinations of relations like `anydo`, so
we would like to speed these up. We would also like to make more calls
terminate. We are interested in improving `simpo`, possibly implementing
Knuth-Bendix Completion(Dick, 1991) relationally. We have done some
preliminary work making a relation to replace `do` that automatic
differentiates forwards and backwards.

##  5\. Related Work

Expresso (Schünemann, 2017) is a computer algebra system written in Clojure
using the miniKanren-inspired core.logic library. expresso’s original intent
was to be relational, but the author made it non-relational to include more
advanced features.(Schünemann, 2020) Like _dxo_ , it includes algebraic
simplification, differentiation, and evaluation. Beyond _dxo_ , it includes
rewriting in normal form and expressions like $sin$.

The Reduce-Algebraic-Expressions system in Prolog (Jasoria, 2019) is similar
to `simpo`, using certain simplification identities. simplifies expressions
like $((x+x)/x)*(y+y-y)$ $\Rightarrow$ $2*y$. It can make simplifications like
$x+x$ $\Rightarrow$ $2*x$ which `simpo` cannot since `simpo` currently only
includes simplification rules involving $0$ and $1$. The Reduce-Algebraic-
Expressions system is not relational.

##  6\. Conclusion

_dxo_ applies relational programming to algebra and differentiation. It can
differentiate, integrate, simplify, complicate, evaluate, create, and reorder.
_dxo_ can concisely represent non-trivial math problems and find solutions.
_dxo_ is a foundation for future exploration of relational programming in
algebra.

###### Acknowledgements.

We are grateful for the work of all the relational programmers whose work we
have built upon. Our work would not have existed without the efforts of Dan
Friedman and Oleg Kiselyov. In addition, we would like to thank Maik
Schünemann for explaining his work. We would like to thank Brandon T. Willard
and Alan T. Sherman for comments on a draft of this paper and sharing ideas
with us. Finally, we thank the anonymous reviewers for their many helpful
comments.

## References

  * (1)
  * Byrd (2009) William E. Byrd. 2009\.  _Relational Programming in miniKanren: Techniques, Applications, and Implementations_.  Ph.D. Dissertation. Indiana University.
  * Byrd et al. (2017) William E. Byrd, Michael Ballantyne, Gregory Rosenblatt, and Matthew Might. 2017.  A Unified Approach to Solving Seven Programming Problems (Functional Pearl).  _Proceedings of the ACM on Programming Languages_ 1, ICFP (2017), 8.
  * Byrd and Friedman (2006) William E. Byrd and Daniel P. Friedman. 2006.  From Variadic Functions to Variadic Relations: A miniKanren Perspective. In _Proceedings of the 2006 Scheme and Functional Programming Workshop_ _(University of Chicago Technical Report TR-2006-06)_ , Robby Findler (Ed.). 105–117.
  * Byrd et al. (2012) William E. Byrd, Eric Holk, and Daniel P. Friedman. 2012\.  miniKanren, Live and Untagged: Quine Generation via Relational Interpreters (Programming Pearl). In _Proceedings of the 2012 Annual Workshop on Scheme and Functional Programming_ _(Scheme ’12)_. ACM, New York, NY, USA, 8–29.  https://doi.org/10.1145/2661103.2661105
  * Dick (1991) A. J. J. Dick. 1991\.  An Introduction to Knuth–Bendix Completion.  _Comput. J._ 34, 1 (01 1991), 2–15.  https://doi.org/10.1093/comjnl/34.1.2 arXiv:https://academic.oup.com/comjnl/article-pdf/34/1/2/1181698/340002.pdf
  * Friedman et al. (2018) Daniel P. Friedman, William E. Byrd, Oleg Kiselyov, and Jason Hemann. 2018.  _The Reasoned Schemer_ (2nd ed.).  The MIT Press, Cambridge, MA, USA.
  * Jasoria (2019) Shruti Jasoria. 2019\.  Reduce-Algebraic-Expressions.  https://github.com/SJasoria/Reduce-Algebraic-Expressions.  Accessed: 2020-07-14.
  * Kiselyov et al. (2008) Oleg Kiselyov, William E. Byrd, Daniel P. Friedman, and Chung-chieh Shan. 2008.  Pure, declarative, and constructive arithmetic relations (declarative pearl). In _Proceedings of the 9th International Symposium on Functional and Logic Programming_ _(LNCS)_ , Jacques Garrigue and Manuel Hermenegildo (Eds.), Vol. 4989. Springer, 64–80.
  * Schünemann (2017) Maik Schünemann. 2017\.  expresso.  https://github.com/clojure-numerics/expresso.  Accessed: 2020-07-14.
  * Schünemann (2020) Maik Schünemann. 2020\.  private correspondence.

##  Appendix A Full Implementation of _dxo_

    (require "faster-miniKanren/mk.rkt")
    (require "faster-miniKanren/numbers.rkt")

    ;defines ^ as expt
    (define ^ (lambda (a b) (expt a b)))

    ;defines ZERO and ONE
    (define ZERO ‘(num ,(build-num 0)))
    (define ONE ‘(num ,(build-num 1)))

    ;exponent: (^ a b)=c
    (define expo
      (lambda (a b c)
        (fresh (bm1 rec)
          (conde
            ((== (build-num 0) b) (== (build-num 1) c))
            ((=/= (build-num 0) b)
             (pluso bm1 (build-num 1) b)
             (expo a bm1 rec)
             (*o a rec c))))))

    ;atom?
    (define atom?
      (lambda (expr)
        (cond
          ((list? expr) #f)
          ((null? expr) #f)
          (else #t))))

    ;atom, null, or list
    (define typeo
      (lambda (expr answer)
        (fresh (a b)
          (conde
            ((== ‘(,a . ,b) expr) (== ’list answer) (=/= ’num a) (=/= ’var a))
            ((== ‘() expr) (== ’null answer))
            ((== ‘(num ,a) expr) (== ’atom answer))
            ((== ‘(var ,a) expr) (== ’atom answer))
            ((== ’+ expr) (== ’+or* answer))
            ((== ’* expr) (== ’+or* answer))))))

    ;miniKanren number
    (define numo
      (lambda (n)
        (fresh (b rest)
          (conde
            ((== ’() n))
            ((== ‘(,b . ,rest) n)
             (conde
               ((== 1 b))
               ((== 0 b)))
             (numo rest))))))

    ;empty env
    (define empty-env ‘())

    ;ext-env
    (define ext-env
      (lambda (x v env)
        (cons ‘(,x . ,v) env)))

    ;lookupo
    (define lookupo
      (lambda (x env v)
        (fresh (env* y w)
          (conde
            ((== ‘((,x . ,v) . ,env*) env))
            ((== ‘((,y . ,w) . ,env*) env) (=/= y x) (lookupo x env* v))))))

    ;unbuild-numinner to every element and if list, then to list
    (define unbuild-numhelper
      (lambda (expr)
        (cond
          ((null? expr) ’())
          ((list? (car expr)) (cons (unbuild-numinner (car expr)) (unbuild-numhelper (cdr expr))))
          (else (cons (car expr) (unbuild-numhelper (cdr expr)))))))

    ;calls unbuld-numinner for every answer in miniKanren
    (define unbuild-num
      (lambda (expr)
        (cond
          ((null? expr) ’())
          (else (cons (unbuild-numinner (car expr)) (unbuild-num (cdr expr)))))))

    ;undoes build-num by calling unbinary
    (define unbuild-numinner
      (lambda (expr)
        (match expr
          [‘() ‘()]
          [‘(num ,b) ‘(num ,(unbinary b 1))]
          [‘(var ,b) ‘(var ,b)]
          [‘(+ ,e . ,e*) (unbuild-numhelper ‘(,e . ,e*))]
          [‘(* ,e . ,e*) (unbuild-numhelper ‘(,e . ,e*))]
          [‘(^ ,e . ,e*) (unbuild-numhelper ‘(,e . ,e*))])))

    ;helper for unbuild-numinner that goes from binary to base 10
    (define unbinary
      (lambda (expr n)
        (cond
          ((null? expr) 0)
          ((atom? expr) expr)
          ((equal? (car expr) 1) (+ n (unbinary (cdr expr) (* 2 n))))
          ((equal? (car expr) 0) (unbinary (cdr expr) (* 2 n))))))

    ;l contains item
    (define membero
      (lambda (item l)
        (fresh (a rest)
          (conde
            ((== ‘(,item . ,rest) l))
            ((== ‘(,a . ,rest) l) (=/= item a) (membero item rest))))))

    ;l does not contain item
    (define notmembero
      (lambda (item l)
        (fresh (a rest)
          (conde
            ((== ’() l))
            ((== ‘(,a . ,rest) l)
             (=/= a item)
             (notmembero item rest))))))

    ;removes item from l
    (define removeo
      (lambda (item contain removed)
        (fresh (rest rest2 c)
          (conde
            ((== ‘(,item . ,rest) contain)
             (== rest removed))

            ((== ‘(,c . ,rest) contain)
             (=/= item c)
             (== ‘(,c . ,rest2) removed)
             (removeo item rest rest2))))))

    "SIMPLIFY";----------------------------------------------------------

    (define simpo
      (lambda (comp simp)
        (fresh ()
          (conde
            ((fresh (n)
               (== ‘(num ,n) comp)
               (== comp simp)))
            ((fresh (v)
               (== ‘(var ,v) comp)
               (== comp simp)))

            ((fresh (e1 e2 s1 s2)
               (== ‘(^ ,e1 ,e2) comp)
               (conde
                 ((== ONE s1) (== ONE simp))
                 ((== ZERO s1) (=/= ZERO s2) (== ZERO simp))
                 ((=/= ZERO s1) (== ZERO s2) (=/= ONE s1) (== ONE simp))
                 ((=/= ZERO s1) (=/= ONE s1) (== ONE s2) (== s1 simp))
                 ((== ‘(^ ,s1 ,s2) simp)
                  (=/= ONE s1)
                  (=/= ONE s2)
                  (=/= ZERO s1)
                  (=/= ZERO s2)))
               (simpo e1 s1)
               (simpo e2 s2)))

            ((fresh (e e* s temp t* n v)
               (== ‘(* ,e . ,e*) comp)
               (conde
                 ((== ’() e*) (simpo e simp))
                 ((== ZERO s)(=/= ’() e*)(== ZERO simp))
                 ((== ONE s)(=/= ’() e*) (simpo ‘(* . ,e*) simp))
                 ((=/= ONE s)
                  (=/= ZERO s)
                  (=/= ’() e*)
                  (conde
                    ((== ZERO temp) (== ZERO simp))
                    ((== ONE temp) (== s simp))
                    ((== ‘(^ . ,n) temp) (== ‘(* ,s ,temp) simp))
                    ((== ‘(+ . ,n) temp) (== ‘(* ,s ,temp) simp))
                    ((== ‘(num ,n) temp) (=/= ZERO temp)
                     (=/= ONE temp) (== ‘(* ,s ,temp) simp))
                    ((== ‘(var ,v) temp) (== ‘(* ,s ,temp) simp))
                    ((==‘(* . ,t*) temp) (== ‘(* ,s . ,t*) simp)))
                  (simpo ‘(* . ,e*) temp)))
               (simpo e s)))

            ((fresh (e e* s temp t* n v)
               (== ‘(+ ,e . ,e*) comp)
               (conde
                 ((== ’() e*) (simpo e simp))
                 ((== ZERO s)(=/= ’() e*)(simpo ‘(+ . ,e*) simp))
                 ((=/= ZERO s)
                  (=/= ’() e*)
                  (conde
                    ((== ZERO temp) (== s simp))
                    ((== ‘(^ . ,n) temp) (== ‘(+ ,s ,temp) simp))
                    ((== ‘(* . ,n) temp) (== ‘(+ ,s ,temp) simp))
                    ((== ‘(num ,n) temp) (=/= ZERO temp) (== ‘(+ ,s ,temp) simp))
                    ((== ‘(var ,v) temp) (== ‘(+ ,s ,temp) simp))
                    ((== ‘(+ . ,t*) temp) (== ‘(+ ,s . ,t*) simp)))
                  (simpo ‘(+ . ,e*) temp)))
               (simpo e s)))))))

    "DERIVATIVE";----------------------------------------------------------

    ;takes derivative
    (define do
      (lambda (x expr deriv)
        (fresh ()
          (symbolo x)
          (conde
            ((fresh (d* e* a b c d)
               (== expr ‘(+ . ,e*)) (== e* ‘(,a . ,b))
               (== deriv ‘(+ . ,d*)) (== d* ‘(,c . ,d))
               (samelengtho e* d*)
               (map-do-o x e* d*)))

            ((fresh ()
               (== expr ‘(^ (var ,x) (num ,(build-num 0))))
               (== deriv ZERO)))

            ((fresh (l e)
               (== expr ‘(* . ,l))
               (letrec ((multruleo
                         (lambda (l dd)
                           (fresh (e e* d d* a b)
                             (conde
                               [(== l ‘(,e))
                                (do x e dd)]
                               [(== l ‘(,e . ,e*))
                                (== e* ‘(,a . ,b))
                                (== dd ‘(+ (* ,d . ,e*) (* ,e ,d*)))
                                (do x e d)
                                (multruleo e* d*)])))))
                 (multruleo l deriv))))

            ((fresh (int intm1)
               (== expr  ‘(^ (var ,x) (num ,int)))
               (== deriv ‘(* (^ (var ,x) (num ,intm1)) (num ,int)))
               (minuso int (build-num 1) intm1)))

            ((fresh (int1 int2)
               (== expr ‘(^ (num ,int1) (num ,int2)))
               (== deriv ZERO)
               (conde
                 ((poso int1))
                 ((== ZERO int1)(poso int2)))))

            ((fresh ()
               (== expr ‘(var ,x))
               (== deriv ONE)))

            ((fresh (int)
               (== expr ‘(num ,int))
               (== deriv ZERO)))))))

    ;maps do relation
    (define map-do-o
      (lambda (x expr* output)
        (fresh (e* e out out*)
          (conde
            [(== expr* ’()) (== output ’())]
            [(== expr* ‘(,e . ,e*))
             (== output ‘(,out . ,out*))
             (do x e out)
             (map-do-o x e* out*)]))))

    "EVALUATE";------------------------------------------------------------

    ;evaluater
    (define evalo
      (lambda (env expr value)
        (fresh (m x c a b rest evc evrest eva evb)
          (conde
            ((== ‘(var ,x) expr) (lookupo x env value))
            ((== ‘(num ,m) expr) (numo m) (== m value))
            ((== ‘(+ ,c . ,rest) expr)
             (conde
               ((== ’() rest) (evalo env c value))
               ((=/= ’() rest)
                (evalo env c evc)
                (evalo env ‘(+ . ,rest) evrest)
                (pluso evc evrest value))))
            ((== ‘(* ,c . ,rest) expr)
             (conde
               ((== ’() rest) (evalo env c value))
               ((=/= ’() rest)
                (evalo env c evc)
                (evalo env ‘(* . ,rest) evrest)
                (*o evc evrest value))))
            ((== ‘(^ ,a ,b) expr)
             (evalo env a eva)
             (evalo env b evb)
             (expo eva evb value))))))

    "REORDER";______________________________________________________________

    ;another option instead of using reordero is to always enter expressions in the same right order
    ;reorders expression deeply, reordering any + and * expressions
    (define reordero
      (lambda (e1 e2)
        (fresh ()
          (conde
            ((== e1 e2) (typeo e1 ’atom))
            ((fresh (o e1* e2*)
               (== ‘(,o . ,e1*) e1)
               (== ‘(,o . ,e2*) e2)
               (typeo o ’+or*)
               (reorderitemso e1* e2*)))
            ((fresh (a1 b1 a2 b2)
               (== ‘(^ ,a1 ,b1) e1)
               (== ‘(^ ,a2 ,b2) e2)
               (reordero a1 a2)
               (reordero b1 b2)))))))

    ;permutes a list by calling reorderinnero, and calls reordero on the items in the list deeply
    (define reorderitemso
      (lambda (e1* e2*)
        (fresh ()
          (samelengtho e1* e2*)
          (reorderinnero e1* e2*))))

    ;permutes and calls reordero on the items, helper for reorderitemso
    (define reorderinnero
      (lambda (e1* e2*)
        (fresh (c1 rc1 rest1 rest2)
          (conde
            ((== ’() e1*)(== ’() e2*))
            ((== ‘(,c1 . ,rest1) e1*)
             (removeo rc1 e2* rest2)
             (reorderinnero rest1 rest2)
             (reordero c1 rc1))))))

    "ANYDO";__________________________________________________________________________________

    (define anydo
      (lambda (expr deriv x)
        (fresh (ecomp dcomp esimp dsimp dcorder)
          (project (expr deriv)
            (if (var? expr)
                (fresh ()
                  (simpo deriv dsimp)
                  (simpo dcorder dsimp)
                  (reordero dcomp dcorder)
                  (do x ecomp dcomp)
                  (simpo ecomp esimp)
                  (simpo expr esimp))

                (fresh ()
                  (simpo expr esimp)
                  (simpo ecomp esimp)
                  (do x ecomp dcomp)
                  (reordero dcomp dcorder)
                  (simpo dcorder dsimp)
                  (simpo deriv dsimp)))))))

    (define doitallevalo
      (lambda (ieval inte deriv deval x env)
        (fresh (icomp dcomp isimp dsimp dcorder)
          (evalo env inte ieval)
          (evalo env deriv deval)
          (do x icomp dcomp)
          (reordero dcomp dcorder)
          (simpo deriv dsimp)
          (simpo dcorder dsimp)
          (simpo inte isimp)
          (simpo icomp isimp))))


# Quantifying Gerrymandering With Simulated Annealing

Stuart Wayland Department of Computer Science, University of California, Santa
Barbara. Research for the completion of the Distinguish in the Major Program.
email<EMAIL_ADDRESS>Faculty Advisor: Professor Eric Vigoda

(August 2022)

###### Abstract

Gerrymandering is the perversion of an election based on manipulation of
voting district boundaries, and has been a historically important yet
difficult task to analytically prove. We propose a Markov Chain Monte Carlo
with Simulated Annealing as a solution for measuring the extent to which a
districting plan is unfair. We put forth promising results in the successful
application of redistricting chains for the state of Texas, using an
implementation of a redistricting Markov Chain with Simulated Annealing to
produce accelerated results. This demonstrates strong evidence that Simulated
Annealing is effective in quickly generating representative voting
distributions for large elections, and furthermore capable of indicating
unfair bias in enacted districting plans.

##  1 Introduction

Legislative redistricting is a necessary and hugely problematic component of
democratic processes. It not only plays an integral role in valuing citizens’
votes, but also has strong implications in the government’s representation of
it’s people. Unfortunately, the choice of a redistricting, and thus the way we
value certain citizens’ votes, may be exploited for partisan gain. This notion
of manipulating district boundaries to favor an outcome, known as
gerrymandering, has become notoriously intertwined with democratic electoral
processes.

In order to comment on the existence of gerrymandering in districting plans,
we must first understand the procedure in which districtings are assigned. For
US Congressional districts, the US Constitution mandates that redistricting
occur every 10 years following the census count and reapportionment of seats
in the US House of Representatives. Furthermore, the Federal Government
requires that each state’s congressional districting plans have equal
population, are contiguous and compact, and comply with the Voting Rights Act
[VRA65].

While these requirements are set forth by the Federal Government, each state
has an elected board of partisan legislators that determines the enacted plan
within these requirements. Should a plan be suspected of being gerrymandered,
it is presented to the respective state’s supreme court who pass judgement on
the fairness of the plan. The reality of judging a plan’s fairness, however,
has been a highly contested and historically difficult task. In past trials,
the efficiency gap has served as the primary indicator for determining a
plan’s fairness. However, as detailed in Chambers et al. [CMS17], this metric
can yield inaccuracies. We will employ a more robust method of evaluating
districting fairness that uses the region’s geo-politics to quantify the
extent to which a districting is gerrymandered. We propose a Markov Chain
Monte Carlo (MCMC) solution to determine districting fairness, as pioneered by
Fifield et al. [FHIT20], Mattingly’s Quantifying Gerrymandering Group at Duke
University [HKL+20], the MGGG at Tuft’s University [DD19], and CMU’s Wesley
Pedgen [PPY17].

An MCMC solution provides a well-principled method for handling the vast
number of potential districtings that adhere to the Federal governments
restrictions on compactness and population equality– we will call these
candidates viable plans. Barring access to corporation sized computational
facilities, enumerating the set of viable plans is computationally
ineffective. Instead, we use an MCMC to sample from the space of viable plans
via a random walk that prefers compact and equal population districtings. We
then apply voting data to the collection of generated plans to produce a
distribution of feasible election results.

This distribution provides an associated probability to the enacted plan that
represents the likelihood the plan should be chosen. In the case that the
enacted plan is gerrymandered, we would expect this plan to have a
significantly low probability of being selected. For example, Professor Wesley
Pegden showed that the 2010 congressional districting of Pennsylvania gave
more seats to Republicans than 99.99% of his randomly generated viable
districting plans, meaning that this plan had a roughly 0.0001 probability of
being fairly selected [PPY17].

Since the emergence of this field in 2018, the previously mentioned MCMC
solution has been presented in court to analyze the presence of
gerrymandering. Two of the aforementioned groups provided expert witness
testimony. Professor Wesley Pegden of CMU testified in Pennsylvania in 2018,
successfully ruling the 2010 Pennsylvania Congressional districting plan as
gerrymandered [Duf18]. Similarly, Professor Johnathan Mattingly of Duke
University provided testimony in North Carolina in 2019, advocating that the
2012 and 2016 North Carolina congressional districting plans were
gerrymandered [Mat21].

We hope to replicate and enhance the experiments and analytical process these
groups employed, expanding this application to the state of Texas. We put
forth promising results in the successful application of redistricting chains
on larger geographic areas while demonstrating a necessity for speed
improvements. We proceed to use a modification of Dr. Benjamin Fifield’s
implementation [KMF+22] of a redistricting Markov Chain to produce accelerated
results using Simulated Annealing, showing strong evidence that this practical
method is effective in quickly generating representative voting distributions.

##  2 Background

This section will give a brief overview of the US Census Bureau’s geographic
division of the United States and explain how the redistricting process can be
framed as a graph cut problem, providing context for an MCMC implementation.
We hope to: 1) provide an understanding of how random districting samples are
generated, and 2) illustrate opportunities to improve districting via a
process called Simulated Annealing.

In the scope of an election process, the US Census Bureau divides each state
into counties, census tracts, and voting tabulation districts (VTDs). VTDs are
the finest geographic division used for aggregating votes, and will be used as
the geographic units in the implementation of our MCMC.

###  2.1 Redistricting as a Graph Cut Problem

Given a shapefile $S$ that divides a state into geographic units, we consider
the dual of the shapefile $G=(V,E)$. The dual $G$ is computed by associating a
vertex with each open face in $S$, and an edge between two vertices if their
corresponding open faces share a border. Also, we have a function
$pop(v):V\to\mathbb{Z}$, where for each $v\in V$,  $pop(v)$ returns the
population of the geographic unit that $v$ represents.

For example, let us consider for simplicity the county shapefile of the state
of Iowa and its corresponding dual graph, as shown in Figure 1. In the right
image, Iowa’s 99 counties are represented by vertices and their borders
represented with edges. Moreover, we can extract a county’s population using
$pop(\cdot)$ on the corresponding vertex.

Figure 1: Iowa County Shape and Corresponding Dual

In the resulting planar graph, we can represent $n$ congressional districts
with an n-coloring. Returning to our Iowa example, we consider representing
Iowa’s 4 congressional districts by coloring each node of the dual graph with
one of four colors, as seen in Figure 2.

Figure 2: Iowa County Dual with Districting Plan

Formally, a districting plan $\sigma:V\to\mathbb{N}$ is an n-coloring of a
graph $G=(V,E)$. For a districting plan $\sigma$, we may also consider
removing all non-Monochromatic edges (edges who have endpoints with different
colors) from G, producing a sub-graph that isolates the state’s districts.
This associates a unique cut edge set $S_{\sigma}$, to each districting
$\sigma$, where $S_{\sigma}\subseteq E$ are the cuts required to separate
colored components in the coloring $\sigma$.

Federal regulations require that districting plans, and therefore graph
colorings, be contiguous. Thus, we can define a valid plan as follows: a
districting plan $\sigma$ is valid if there exists a cut edge set
$S_{\sigma}\subseteq E$ such that the sub graph, $(V,E\setminus S_{\sigma})$
has exactly $n$ connected components $C_{\sigma}^{1}$, …, $C_{\sigma}^{n}$.
Let $\Omega$ be the set of all valid plans.

These connected components $C_{\sigma}^{1},...,C_{\sigma}^{n}$ are the
districts laid out by the plan $\sigma$, and will be referred to as such.
Thus, we can calculate the population of a district by aggregating the
population of each vertex in the connected component. That is, given a
district $C_{\sigma}^{i}\subseteq V$, we define the population of
$C_{\sigma}^{i}$ as $pop(C_{\sigma}^{i})=\sum_{v\in C_{\sigma}^{i}}pop(v).$ It
then follows that $pop(V)$ is equal to the total population of the state.

Now, we can define two important metrics on a districting plan $\sigma$ and
it’s cut edge set $S_{\sigma}$. First, the population equality metric:

$pop\\_eq(\sigma)=\sum_{i=1}^{n}\left|pop(C_{\sigma}^{i})-\frac{pop(V)}{n}\right|$

And next, the compactness metric, which is determined by the ratio of cut
edges to total edges in the dual graph G=(V,E):

$comp(\sigma)=\frac{|S_{\sigma}|}{|E|}.$

We hope to favor a smaller subset of $\Omega$ that obey Federal regulations
for compactness and population difference. These will be referred to as viable
plans, and will be targeted by using the score functions $comp()$ and
$pop\\_eq$.

We will generate a viable plan by performing a random walk on $\Omega$ as
follows: suppose we have a valid districting plan $\sigma\in\Omega$, defined
by the edge set $S_{\sigma}$ for the dual graph $G=(V,E)$. We begin by
selecting a random vertex $v\in V$ such that $N(v)\cap
S_{\sigma}\not=\emptyset$. We then select a random edge $e\in N(v)\cap
S_{\sigma}$. That is, $v$ is a vertex lying on an internal boundary of a
district. We then recolor $v$ to the neighboring district which $v$ is
connected to through $e\in N(v)\cap S_{\sigma}$. Finally, we check that the
resulting districting plan from this transition, call it $\sigma^{\prime}$,
belongs to the space of valid districtings, $\Omega$. If $\sigma^{\prime}$ is
not a valid districting, we repeat the process on $\sigma$.

Once we have produced a new valid districting $\sigma^{\prime}$, we apply a
Metropolis Filter. This will leverage the metrics discussed previously to
favor plans with more desirable compactness and population equality
properties. This filter sets the probability of transitioning to a new plan
equal to the weight of the new plan over the weight of the old. More
precisely, for a weight function $w:\Omega\to\mathbb{R}$, we set the
transition probability,

$p=\min\left\\{\frac{w(\sigma^{\prime})}{w(\sigma)},1\right\\}.$

The Metropolis Filter guarantees that our sampling probability distribution is
proportional to our weight function $w$, allowing us to control the
districting plans we wish to sample. With this in mind, we set our weight
function as

$w(\sigma)=e^{-(\beta_{1}pop\\_eq(\sigma)+\beta_{2}comp(\sigma))},$

where $\beta_{1},\beta_{2}$ are constants. We transition to the plan
$\sigma^{\prime}$ with our above probability $p$, setting
$\sigma=\sigma^{\prime}$, or otherwise leave $\sigma$ unchanged. We then
repeat the process until we have arrived at a randomly generated districting
plan. The number of steps required to achieve a random sample is known as the
chain’s mixing time [JS96], and is a theoretical topic altogether unknown in
the scope of redistricting chains. Our goal is to sample from a distribution
proportional to $w$ with very large $\beta_{1},\beta_{2}$ values. However,
these large $\beta_{1},\beta_{2}$ values severely constrict the mixing of a
chain over the space $\Omega$, whereas small $\beta_{1},\beta_{2}$ allow
greater movement over space and therefore faster apparent mixing. For this
reason, MCMC chains often employ methods to artificially speed up mixing in
practice, even if these methods do not ensure any theoretical difference in
the mixing time.

###  2.2 Simulated Annealing

One such practical method is known as simulated annealing, which dynamically
adjusts the $\beta_{1},\beta_{2}$ values as the chain is mixing. The
$\beta_{1},\beta_{2}$ values control the influence of population equality and
the compactness metric respectively. This allows us to manage how severely we
wish to enforce our constraint metrics. For example, consider
$\beta_{1},\beta_{2}=0$, meaning for all $\sigma\in\Omega$, $w(\sigma)=1$.
This removes all influence of our weight function, resulting in uniformly
sampling from our space $\Omega$ of valid plans. Conversely, if we consider
extremely large $\beta_{1},\beta_{2}$, then our weight function approaches
zero for all plans without a zero population equality (all districts have
exactly the same population) and compactness score. This results in only
sampling extremely compact, equal population districting plans. In the context
of a normal Markov Chain, we choose our $\beta_{1},\beta_{2}$ values so that
our sampled plans fall within federal laws of compactness and population
equality. We will refer to the $\beta_{1},\beta_{2}$ values that produce such
plans as the target values.

To provide context to the name simulated annealing, we may think of
$\beta_{1},\beta_{2}$ values as inverse temperature. Recall that annealing is
the process of heating a metal before allowing it to cool slowly in a desired
shape. Similarly, simulated annealing refers to beginning a Markov chain with
$\beta_{1},\beta_{2}=0$. This may be thought of as heating the chain, and is
often referred to as the "hot steps", in which the chain may move freely
across the entire space of valid districtings. In order to slowly "cool" the
chain, the $\beta_{1},\beta_{2}$ values are increased at a rate called the
cooling schedule until they reach their target values. The final steps of the
chain at the target $\beta_{1},\beta_{2}$ values are known as the "cold
steps", and ensure we are sampling from the desired distribution guaranteed
from the Metropolis Filter.

The benefit of simulated annealing is its ability to increase the speed of
plan generation. This improvement comes from decreasing the total number of
steps necessary to produce a near-random plan. Since the majority of the chain
is run with small $\beta_{1},\beta_{2}$ values, the probability of accepting a
proposed plan in the chain is significantly higher than its pure Markov chain
counterpart. This means we may achieve the same number of accepted steps, and
so in practice the same amount of "randomness", with far fewer iterations of
the chain.

##  3 Results

We generated randomly sampled districting ensembles for the states of North
Carolina and Texas, using North Carolina as a small scale validation study for
our simulated annealing algorithm. We compared our results in North Carolina
to the findings of Herschlag et al. [HKL+20].

We analyze the sets of generated plans, which we will refer to as ensembles,
by processing election outcomes for each plan in the ensemble. These election
outcomes produce a distribution that we compare to an enacted plan. From the
selection of our weight function, we expect to see a near-normal distribution
of election outcomes. As such, we fit a density curve to our election outcome
data to provide an associated probability to our enacted plan, which we assume
to be the probability this outcome should occur fairly.

###  3.1 North Carolina

We begin by sampling an ensemble of congressional districting plans for the
2016 North Carolina US House of Representatives election. We will be using the
2010 North Carolina VTD shapefile, and will be comparing to the North Carolina
congressional districting enacted in 2010.

We first replicated the study done by Herschlag et al. [HKL+20], producing
similar results before testing our simulated annealing algorithm. We generated
50 plans with the Redist flip algorithm to validate our preprocessing and
analytic procedures. We compared our distribution of 50 plans to the 66,544
plans generated by Herschlag et al. [HKL+20] We then fitted density curves to
both distributions, as seen in Figure 3.

Figure 3: Distribution of number of Democratic seats won among the 13 seats in
North Carolina. We display the results of Duke University (left) and our
Redist Flip algorithm (right). In both figures, the red line represents the
outcome of the enacted plan.

Despite only generating 50 plans, we see a similarity to the results in
Herschlag et al. [HKL+20], who showed plans giving 3 Democratic seats in this
election had a probability of 0.008 of being fairly selected. Unsurprisingly,
we did not produce a single plan yielding 3 Democratic seats. Regardless of
our inadequate number of generated plans, we did validate that our
preprocessing of voting/geographic data and our election analysis both align
with previous work. We fitted a normal density curve to our plan distribution,
calculating a z-score for 3 Democratic seats. We found that a plan yielding a
3 Democratic seat outcome had an approximate probability of 0.0021 of being
selected.

Figure 4: Distribution of number of Democratic seats won among the 13 seats in
North Carolina. We display the results of Duke University (left) and our
simulated annealing algorithm (right). In both figures, the red line
represents the outcome of the enacted plan.

After validating our data processing and analysis methods, we generated 500
plans for North Carolina using simulated annealing. As shown in Figure 4, we
obtained a similar distribution as Duke, with the vast majority of sampled
plans yielding 4, 5, and 6 Democratic seats in Congress, compared to the 3
seats awarded in the 2016 election. Like our flip algorithm, not a single plan
generated with simulated annealing produces an outcome of 3 Democratic seats
won. Once again, we fitted a normal density curve to our generated
distribution to calculate the probabilities of producing an election outcome
with 3 Democratic seats won. Our calculated probabilities compared to
Mattingly’s Duke group[HKL+20] can be seen in Table 1. These results show
promise in expanding simulated annealing to a far larger application, namely
plan generation for the state of Texas.

Democratic Seats Won | 3 | 4 | 5 | 6
---|---|---|---|---
Redistricting Chain | 0.0021 | 0.21 | 0.46 | 0.38
Simulated Annealing | 0.0013 | .33 | .61 | .06
Duke Group | 0.008 | .29 | .52 | .11
Table 1: Approximate probabilities of 2016 NC USH Election Outcomes from
Redistricting chain, Simulated Annealing Distribution and Duke’s Quantifying
Gerrymandering Group[HKL+20]

###  3.2 Texas

Like in North Carolina, we produce a voting distribution of an ensemble of
generated districting plans to analyze outcomes of the 2020 US House of
Representatives vote. We will be comparing our generated plans to the 2020
Texas Congressional districting, and again will be producing a distribution
from the Redist flip algorithm as well as a distribution generated from
simulated annealing. We generated 50 plans using the Redist algorithm and 500
plans using simulated annealing. These two distributions can be seen in Figure
5.

Figure 5: Distribution of number of Democratic seats won in 2020 US House of
Representatives election in Texas. Shown left is distribution of 50 plans
generated from Redist Flip algorithm, and the right shows distribution of 500
plans generated with simulated annealing. The vertical red line represents the
real outcome from the enacted Congressional Districting plan.

Despite the 2020 Texas Congressional Districting plan yielding 13 of 36 seats
to Democrats, our Redist flip algorithm, shown on the left of Figure 5, did
not produce a single plan with 13 or even 14 Democratic seats won. The
majority of plans generated gave between 16 and 18 seats to Democrats.

Similarly, our simulated annealing algorithm also did not produce a single
plan with the same 13 Democratic seat outcome. Like the flip algorithm, the
majority of the 500 generated plans gave between 16 and 18 seats to Democrats.
We fitted normal density curves to both distributions, calculating the
probability of each election outcome from our simulated annealing
distribution, as shown in Table 2.

Democratic Seats Won | 13 | 14 | 15 | 16 | 17 | 18 | 19
---|---|---|---|---|---|---|---
Probability | 0.0002 | 0.01 | .11 | .29 | 0.44 | .13 | .03
Table 2: Approximate probabilities of 2020 TX USH Election Outcomes from
Simulated Annealing Distribution

##  4 Methodology

In this section, we describe the proposed methodology for running the Redist
package’s flip algorithm in the states of North Carolina, Iowa, and Texas. We
will provide the specifications of the Markov Chain algorithm used to generate
the voting distributions, the preprocessing methods of external shape files
and voting data sets, and the precise simulation inputs and reasoning behind
the parameter choices. Considering the state-specific data sources, these
preprocessing and parameter choices will vary based on the size of the state
and the manner in which voting and population data is recorded. Furthermore,
we will detail the simulated annealing algorithm used to improve the speed of
plan generation, reporting the specifications yielding the results provided
above.

###  4.1 The Redist Flip Algorithm

The Redist flip algorithm is an expansion of the Swendsen-Wang Algorithm
[BZ05], a simple Markov Chain described in Section 2. Now, let us set about
describing the chain used to generate the above results.

We will then set forth describing the transition on step $t$ of the flip
chain, wherein we obtain $\sigma_{t+1}$, defined by edge set
$S_{\sigma_{t+1}}$ from $\sigma_{t}$, defined by edge set $S_{\sigma_{t}}$.

Seeing how $S_{\sigma}$ refers to the cut, or non-monochromatic edges in the
plan $\sigma$, let us define $\overline{E}(\sigma)$ to be the set of uncut, or
monochromatic edges in the plan $\sigma$. That is,

$\overline{E}(\sigma):=\left\\{e\in E\backslash S_{\sigma}\right\\}.$

The algorithm begins by labelling each edge $e\in\overline{E}(\sigma_{t})$ as
a "flip" edge with probability $\lambda$. Let us refer to this labelling of
edges $e\in\overline{E}(\sigma_{t})$ as $\pi$. We then define

$\overline{E}_{flip}(\sigma_{t},\pi)=\left\\{e\in\overline{E}(\sigma_{t}):\pi(e)=1\right\\}.$

In other words, $\overline{E}_{flip}(\sigma_{t},\pi)$ is the set of "flip"
edges resulting from $\pi$.
$\overline{E}_{flip}(\sigma_{t},\pi)\subseteq\overline{E}(\sigma_{t})$ can
also be thought of as creating a set of connected components lying within
$\overline{E}(\sigma_{t})$. Let us call this set of connected components
resulting from our labeling of flip edges $C_{flip}(\sigma_{t})$. Next, let us
consider the set $C_{\partial}(\sigma_{t})\subseteq C_{flip}(\sigma_{t})$,
which we define as

$C_{\partial}(\sigma_{t}):=\left\\{C\in C_{flip}(\sigma_{t}):\exists v\in
C\text{ s.t. }\exists\text{ }e\in N(v)\text{ where }e\in
S_{\sigma_{t}}\right\\}.$

That is, $C_{\partial}(\sigma_{t})$ is the set of connected components from
the "flip" labeling that lie on the boundary of the plan $\sigma_{t}$. Next,
the algorithm selects a random subset $F\subseteq C_{\partial}(\pi_{t})$ such
that $F$ does not contain any adjacent connected components, and flips each
component within $F$ its respective neighboring coloring. In the case that a
component has multiple neighboring districts, a neighboring district is chosen
at random. Let us call the resulting plan from this flip of components
$\sigma^{\prime}$.

Next, we ascertain the validity of the proposed plan $\sigma^{\prime}$ by
checking each district is contiguous. If $\sigma^{\prime}$ is not a valid
districting, it is rejected and the process is repeated for $\sigma_{t}$.

Finally, this algorithm applies a Metropolis-Hastings filter to determine a
probability of transitioning to $\sigma^{\prime}$ by setting
$\sigma_{t+1}=\sigma^{\prime}$. The algorithm accepts the proposed plan as
$\sigma_{t+1}$ with a probability proportional to its desirability. To do
this, we use the weight function $w(\sigma)$, where the probability of
accepting $\sigma^{\prime}$ as $\sigma_{t+1}$ is

$Pr(\text{accepts $\sigma^{\prime}$ as $\sigma_{t+1}$})=\begin{cases}1&\text{
if }w(\sigma^{\prime})>w(\sigma_{t})\\\
\frac{w(\sigma^{\prime})}{w(\sigma_{t})}&\text{ if }w(\sigma^{\prime})\leq
w(\sigma_{t})\end{cases}.$

As stated above, we define our weight function as

$w(\sigma)=e^{-(\beta_{1}pop\\_eq(\sigma)+\beta_{2}comp(\sigma))}.$

The Flip algorithm continues this process until completing its target number
of accepted steps, where an accepted step is the case where the algorithm
transitions to a proposed plan, setting $\sigma_{t}=\sigma^{\prime}$.

###  4.2 Reproducing a study in North Carolina

In order to investigate the validity of the Redist flip algorithm, we set
about reproducing the results of Herschlag et al. [HKL+20]. For this
experiment, we use the same voting, population, and shapefile data, retrieving
these data sets from the US Census Bureau [Bur]. We will be using the 2010
Voting Tabulation Shape File for North Carolina and the 2016 US House of
Representatives Voting data. We then set about selecting parameters that would
produce the most reliable and efficient simulations.

First, we looked at restricting our plans to an acceptable population range
for each district. In compliance with Supreme Court ruling Karcher v. Daggett
[KD83], the population tolerance of Congressional Districts should be equal
unless the state can provide good reason for the unbalance. Despite this
ruling, running the flip algorithm without a population tolerance is not
feasible, as the acceptance rate becomes negligibly close to 0, making
movement of the chain near impossible. For this reason, we select a population
tolerance of $10\%$, as it provides the smallest population variance while
still complying with Brown v. Thompson, which allows for a $10\%$ variance in
State and local districting elections [AM].

Next, we select a constraint weight for compactness. This parameter, however,
was selected almost entirely by its affect on running time. Since compactness
itself is not precisely defined in political science, we simply aimed to avoid
plans that overtly violate compactness. We found that a compactness score of
.4 results in a manageable running time, while providing plans that satisfy a
loose visual compactness metric. The measuring of compactness itself is a
difficult and problematic field, and is discussed in further detail by Roland
Fryer and Richard Holden [Fry07].

Third, we select the number of simulation steps to run for each sampled plan.
Since there is no theoretical justification of a mixing time for this chain,
we must resign ourselves to using approximation metrics to determine the
mixing of the chain. We used the Gelman-Rubin diagnostic of convergence to
obtain a minimum number of simulation steps required. Using 10 independent
chains, we increased the number of simulation steps until we received a
Gelman-Rubin less than 1.1, as is the indication of a well-mixing chain
[GR92]. With this metric, we decided to run 3000 simulation steps in our
recreation of North Carolina for each sampled plan. We ran this simulation for
6 days to obtain 50 sampled plans, serving as a small ensemble of plans for
North Carolina.

Pop Tolerance | Constraint | Weight | Nsims | Nplans
---|---|---|---|---
0.1 | Edges Removed | 0.4 | 3000 | 50
Table 3: North Carolina Redist Flip Parameters

Finally, we calculated the number of Democratic seats won in the 2016 House of
Representatives Vote, plotting the distribution of Democratic seats won to
compare to results of Herschlag et al. [HKL+20].

###  4.3 Generating a Simulated Annealing Algorithm

Improving the speed of generating each sampled plan is satisfied through
simulated annealing. Working with the 2020 Iowa County level shapefile, we
designed an algorithm to run simulation steps without any constraint
restriction before slowly annealing back to the population and compactness
constraints described above. We choose to use the Iowa County shapefile as a
basis for design because of its simplistic nature, allowing for fast
experimentation on the cooling schedule, or rate at which we reintroduce the
population and compactness constraints.

First, we run 500 steps without constraint weights, which we refer to as hot
steps. These steps contribute to the majority of the mixing of the chain, and
since we are running without constraints, we can complete them quickly with a
100% acceptance rate. When running 1000 hot steps, we see a similar
distribution of sampled plans, allowing us to loosely conclude that 500 hot
steps is sufficiently mixing the chain.

To achieve the same desired distribution as the flip algorithm with
constraints, we need to re-introduce the constraint parameters at a rate that
prevents a near-zero acceptance rate. Reintroducing the constraint too quickly
will prevent the chain from accepting any proposed plan, resulting in a
completely halted, or stagnant chain.

With this goal in mind and after experimentation, we decided to choose a
population tolerance change of 0.005, or .5%. That is, we would run a flip
chain for 10 steps before decreasing the population tolerance by 0.05.
Starting with a population tolerance of 1, we would repeat this procedure
until arriving at our target population tolerance of .1, or 10%. With these
parameter choices, we would require 1800 annealing steps for population alone.

Similarly, starting with a compactness weight of 0, we increase our
compactness constraint weight by 0.01 until reaching our desired weight of
0.4. The result is 400 steps to reintroduce the desired compactness weights.
The choice to anneal population and compactness separately, along with the
cooling schedule were chosen to minimize the number of steps necessary to
return to desirable constraint weights, while avoiding with likely probability
the case of a stagnant chain.

Hot Steps | Anneal Steps | Cold Steps | Pop Tolerance delta | Compact delta
---|---|---|---|---
500 | 10 | 100 | 0.005 | 0.01
Table 4: Simulated Annealing Parameters

Finally, we run the flip algorithm for 100 steps with our desired constraints,
which we will call cold steps. The purpose of these 100 steps is not to
provide any further mixing, but instead to ensure that our resulting plan has
an acceptable compactness and population score.

After configuring our simulated annealing algorithm, we tested it in North
Carolina. We compare a distribution composed of 300 sampled plans from our
annealing algorithm to the distribution produced by Herschlag et al. [HKL+20].

###  4.4 Expanding to Texas

Proceeding in a similar fashion to that of North Carolina, we set about
generating an ensemble for Texas Congressional districting plans. We retrieve
the 2020 Texas VTD shapefile and population data from the Redistricting Data
Hub, which was collected by the US Census Bureau. We also retrieve the 2020 US
House of Representatives election also from the US Census Bureau [Bur].

We then set about preprocessing the 2020 US House data file to align with our
Texas VTD shapefile. We iterate through each county of Texas, matching
precinct names recorded in the voting data to VTD names recorded in our
population data. In the case that multiple precincts fall into the same VTD,
we sum the precinct votes and record it as the vote for the VTD. Once
completing this matching procedure by name, we matched the remaining unmatched
precincts based on geographic data through a Google maps query. After
assigning all precinct voting data to a VTD, we still had 9 out of 9007 VTDs
which had not yet been assigned a vote from the election. These VTD’s were
extremely sparsely populated regions, with all having a population less than
10. The majority (6 of 9) of these problematic VTD’s were located in County 29
and 439, while the other 3 were located in County 37, 41, and 415. We assigned
both the Democratic and republican vote of these VTDs as half the population
of the VTD.

Using the same constraints as in North Carolina, or population tolerance of 10
% and compactness weight of 0.4, we choose the minimum number of simulation
steps required to achieve a Gelman-Rubin score less than 1.1 when running 10
separate chains. We arrive at 5000 necessary steps, and sample 50 plans both
with a compactness constraint and without a compactness constraint.

Pop Tolerance | Constraint | Weight | Nsims | Nplans
---|---|---|---|---
0.1 | Edges Removed | 0.4 | 5000 | 50
Table 5: Texas Redist Flip Parameters

The parameters for the chains creating this ensemble can be found in Table 5.
While more than 50 sampled plans is preferable, the amount of time required to
produce each individual time dictated our 50 plan total.

Taking the same approach for our simulated annealing algorithm in Texas as in
North Carolina, we generated 500 plans using the same parameters shown in  3.

Finally, we plotted the number of Democratic seats won for each ensemble of
plans representing our distributions, comparing to the number of seats won in
the existing plan used in the 2020 US House of Representatives Vote.

##  5 Discussion

We found that the 2020 Texas Congressional Districting gave an unfair
advantage to Republicans in the 2020 US House election. The results shown in
Figure 5 demonstrate an abnormal outcome in this election, which not only
favors Republicans unfairly but is also too statistically improbable to be
unintentional. Our distribution indicates that the enacted plan election
outcome is 3 standard deviations from our generated mean, suggesting that fair
selection of a such a plan is infeasible. These findings not only convey
manipulation within the redistricting system that conventional metrics of
gerrymandering detection have failed to effectively address, but exhibit the
2020 Texas US House congressional election outcome as an extreme outlier of
feasible election outcomes.

We emphasize that we do not see this study as adequately concluding on the
degree to which Texas has an unfair congressional districting. Our validation
study has shown that insufficient samples may lead to distorted probabilities
on election outcomes. Seeing how we were not able to generate a large enough
sample size from which we created our reference distribution, we understand
that our 0.0002 probability the 2020 US House election outcome with a new
districting plan are not representative. The small sample size is likely to
significantly decrease this approximated likelihood, producing an inaccurate
quantification of gerrymandering in Texas.

However, we did produce 550 near-random viable congressional districting plans
for Texas without a single one yielding the result seen in the actual
election. This is clear evidence indicating that a Congressional districting
plan yielding 13 Democratic seats must be intentionally manipulated for the
benefit of the Republican party. While we do not argue that a fair plan would
yield our mean of 17 Democratic seats won, we believe such a drastic
difference should inspire further investigation into the extent to which such
a result is unfair. With study as motivation, we hope to see further
investigation into larger elections with recognition that Simulated Annealing
is a promising method allowing political analysis of larger geographic areas.

## References

  * [AM] Micah Altman and Michael McDonald.  Redistricting Criteria: Equal Population - Public Mapping Project.
  * [Bur] US Census Bureau.  TIGER/Line Shapefiles.  Section: Government.
  * [BZ05] Adrian Barbu and Song-Chun Zhu.  Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities.  IEEE transactions on pattern analysis and machine intelligence, 27(8):1239–1253, August 2005.
  * [CMS17] Christopher P. Chambers, Alan D. Miller, and Joel Sobel.  Flaws in the Efficiency Gap.  Journal of Law & Politics, 33:1, 2017.
  * [DD19] Daryl DeFord and Moon Duchin.  Redistricting Reform in Virginia: Districting Criteria in Context.  Virginia Policy Review, 12(2):22, 2019.
  * [Duf18] Jocelyn Duffy.  Pa. Congressional Map Declared Unconstitutional - News - Carnegie Mellon University, January 2018.
  * [FHIT20] Benjamin Fifield, Michael Higgins, Kosuke Imai, and Alexander Tarr.  Automated Redistricting Simulation Using Markov Chain Monte Carlo.  Journal of Computational and Graphical Statistics, 29(4):715–728, October 2020.
  * [Fry07] Roland G Fryer.  Measuring the Compactness of Political Districting Plans.  National Burea of of Economic Research, page 52, October 2007.
  * [GR92] Andrew Gelman and Donald B. Rubin.  Inference from Iterative Simulation Using Multiple Sequences.  Statistical Science, 7(4):457–472, 1992.  Publisher: Institute of Mathematical Statistics.
  * [HKL+20] Gregory Herschlag, Han Sung Kang, Justin Luo, Christy Vaughn Graves, Sachet Bangia, Robert Ravier, and Jonathan C. Mattingly.  Quantifying Gerrymandering in North Carolina.  Statistics and Public Policy, 7(1):30–38, January 2020.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/2330443X.2020.1796400.
  * [JS96] Mark Jerrum and Alistair Sinclair.  The Markov chain Monte Carlo Method: an approach to approximate counting and integration.  Approximation Algorithms for NP-hard problems, PWS Publishing, 1996\.
  * [KD83] KD.  Karcher v. Daggett, 462 U.S. 725, 1983.
  * [KMF+22] Christopher T. Kenny, Cory McCartan, Ben Fifield, Kosuke Imai, Jun Kawahara, Alexander Tarr, and Michael Higgins.  redist: Simulation Methods for Legislative Redistricting, June 2022.
  * [Mat21] Jonathan C Mattingly.  Expert Report on the North Carolina State Legislature and Congressional Redistricting (Corrected Version), December 2021.
  * [PPY17] Wesley Pegden, Ariel D. Procaccia, and Dingli Yu.  A partisan districting protocol with provably nonpartisan outcomes, October 2017.  arXiv:1710.08781 [cs, math].
  * [VRA65] VRA.  Voting Rights Act, 1965.