## Ray Tracing Algebraic Surfaces

Ray tracing is a natural way of producing computer images. One takes a geometrical ray that connects the pinhole of the camera to a pixel of the camera and find where it hits objects in the scene. You then color the pixel the color of the object it hit.

You can add a great deal of complexity to this by more sophisticated sampling and lighting, multiple bounces, strange surfaces, but that’s it in a nutshell.

A very popular tutorial on this is Ray Tracing in One Weekend https://raytracing.github.io/

There are a couple ways to do the geometrical collision detection part. One is to consider simple shapes like triangles and spheres and find closed form algorithms for the collision point. This is a fast and simple approach and the rough basis of the standard graphics pipeline. Another is to describe shapes via signed distance functions that tell you how far from the object you are and use ray-marching, which is a variant of newton’s method iteratively finding a position on a surface along the ray. ShaderToys very often use this technique.

If you describe your objects using algebraic (polynomial) equations, like $x^2 + y^2 + z^2 - 1$ describes a sphere, there is the possibility of using root finding algorithms, which are readily available. I thought this was kind of neat. Basically the ray hitting the concrete pixel $(x_0, y_0)$ can be parameterized by a univariate polynomial $(x,y,z) = (\lambda x_0, \lambda y_0, \lambda)$ , which can be plugged into the multivariate polynomial $(\lambda x_0)^2 + (\lambda y_0)^2 + \lambda^2 - 1$. This is a univariate polynomial which can be solved for all possible collision points via root finding. We filter for the collisions that are closest and in front of the camera. We can also use partial differentiation of the surface equations to find normal vectors at that point for the purposes of simple directional lighting.

As is, it really isn’t very fast but it’s short and it works.

Three key packages are

using Images
using LinearAlgebra
using TypedPolynomials
using Polynomials

function raytrace(x2,y2,p)
z = Polynomials.Polynomial([0,1])

# The ray parameterized by z through the origin and the point [x2,y2,1]
x3 = [z*x2, z*y2, z]

# get all the roots after substitution into the surface equation
r = roots(p(x=>x3))

# filter to use values of z that are real and in front of the camera
hits = map(real, filter( x -> isreal(x) & (real(x) > 0.0)  , r))

if length(hits) > 0
l = minimum(hits) # closest hit only
x3 = [z(l) for z in x3]
# get normal vector of surface at that point
dp = differentiate(p, x)
normal = normalize([ z(x=> x3)  for z in dp])
# a little directional and ambient shading
return max(0,0.5*dot(normal,normalize([0,1,-1]))) + 0.2
else
return 0 # Ray did not hit surface
end
end

@polyvar x[1:3]

# a sphere of radius 1 with center at (0,0,3)
p = x[1]^2 + x[2]^2 + (x[3] - 3)^2 - 1

box = -1:0.01:1
Gray.([ raytrace(x,y,p) for x=box, y=box ])

Sphere.

@polyvar x[1:3]
R = 2
r = 1

# another way of doing offset
x1 = x .+ [ 0, 0 , -5 ]

# a torus at (0,0,5)
# equation from https://en.wikipedia.org/wiki/Torus
p = (x1[1]^2 + x1[2]^2 + x1[3]^2 + R^2 - r^2)^2 - 4R^2 * (x1[1]^2 + x1[2]^2)

box = -1:0.005:1
img = Gray.([ raytrace(x,y,p) for x=box, y=box ])
save("torus.jpg",img)

Some thoughts on speeding up: Move polynomial manipulations out of the loop. Perhaps partial evaluate with respect to the polynomial? That’d be neat. And of course, parallelize

## A Short Skinny on Relations & the Algebra of Programming

I’ve been reading about the Algebra of Programming lately and lovin’ it. See J.N. Oliveira’s draft text in particular and the links in the references. I’ve started exploring the stuff from this post and more over here: https://github.com/philzook58/rel

## Why and What?

Relations can expand the power of functional programming for the purpose of specification.

The point of a specification is to be able to write down in a very compact and clear way your intent for a program, more clearly and compactly than a full implementation could be written. It therefore makes sense to add to your specification language constructs that are not necessarily executable or efficient for the sake of compactness and clarity. When one needs executability or efficiency, one writes an implementation whose behavior you can connect to the spec via a formal or informal process.

Functional programming, with it’s focus on the abstraction of the mathematical function, is a good trade-off between executability, efficiency, and expressibility. It lies in a reasonable location between the ideas amenable to reasoning by a human mind and the command-driven requirements of the machines.

Functions are a specialization of relations. Relations extend the mathematical notion of functions with constructs like nondeterministic choice, failure and converse. These constructs are not always obviously executable or efficient. However, they greatly extend the abilities of reasoning and the clarity of expression of a specification.

The point-free style of reasoning about functions extends to a point-free style reasoning about relations, which is known as relation algebra. There are rich analogies with databases, category theory, linear algebra, and other topics.

Plus, I think it is very neato for some reason. If anyone ever thinks something is really neato, isn’t it worth giving it a listen?

### A Simple Representation of Relations in Haskell

The simplest description of relations is as a set of tuples. So first let’s talk a bit about the options for sets in Haskell.

There are a couple different reasonable ways to represent sets in Haskell.

• [a] or Vector a
• a -> Bool
• Set a — a tree based Set from the containers package.

These have different performance characteristics and different power. The list [a] is very simple and has specialized pleasant syntax available. The indicator function a -> Bool gives you no ability to produce values of type a, but can easily denote very sophisticated spaces. Set a is a good general purpose data structure with fast lookup. You might also choose to mix and match combinations of these. Interconversion is often possible, but expensive. This is not a complete list of possibilities for sets, for example you may want a representation with a stronger possibility for search.

We can directly use the definition of relations as a set of tuples with the above

type Rel a b = [(a,b)]
type SetRel a b = Set (a,b)
type FunRel a b = (a,b) -> Bool

But we also have the option to “curry” our relation representations, sort of mixing and matching properties of these representations.

type List a b = a -> [b] -- Commonly seen type  in List monad/logic programming
type MapRel a b = Map a (Set b)

You might also choose to package up multiples of these representations, choosing the most appropriate as the situation requires, see for example the relation package, whose type holds both Map a (Set b) and Map b (Set a).

Despite fiendishly poor performance, for simplicity and list comprehension syntax we are going to be using type Rel a b = [(a,b)] for the remainder of the post.

I’m also taking the restriction that we’re working in bounded enumerable spaces for ease. I assume such a requirement can be lifted for many purposes, but finite spaces like these are especially well tamed. The following typeclass and definition is very useful in this case.

type BEnum = (Enum a, Bounded a)
enumAll :: (BEnum a) => [a]
enumAll = [minBound .. maxBound]

#### Functions and Relations

Functions can be thought of as relations with the special property that for each left part of the tuple, there is exactly one right side and every possible left side appears. The relation corresponding to a function $f$ looks like $F = \{(x,y) | x \in X, y \in Y, y = f (x)\}$.

tabulate :: (BEnum a) => (a -> b) -> Rel a b
tabulate f = [(x, f x) | x <- enumAll]

There is a natural and slightly clever lifting of function composition to relations. We now check whether there exists a value that is in the right side of one and the left side of the other.

rcompose :: Eq b => Rel b c -> Rel a b -> Rel a c
rcompose xs ys = [ (a,c)  | (a, b) <- ys, (b', c) <- xs, b' == b]

rid :: (Enum a, Bounded a) => Rel a a
rid = tabulate id

Because of these two operations (and their properties of associativity and absorption), FinRel is a category. We do however need the Eq b restriction to make Rel an instance of the category typeclass, so it does not quite fit the definition of category in base. It is a constrained category.

We can lift the common arrow/categorical combinators up to relations for example.

-- arrow/category combinators
rfan :: Eq a => Rel a b -> Rel a c -> Rel a (b,c)
rfan f g = [ (a, (b,c)) | (a,b) <- f, (a',c) <- g, a == a']

rfst :: BEnum (a,b) => Rel (a,b) a
rfst = tabulate fst

rsnd :: BEnum (a,b) => Rel (a,b) b
rsnd = tabulate snd

rleft :: (Enum a, Bounded a) => Rel a (Either a b)
rleft = tabulate Left

rright :: BEnum b => Rel b (Either a b)
rright = tabulate Right

reither :: Eq a => Rel a c -> Rel b c -> Rel (Either a b) c
reither f g = [(Left a, c) | (a,c) <- f] ++ [(Right b, c) | (b,c) <- g]


With these combinators, you have access to many functions on basic non-recursive algebraic data types. By combining them in a point free style, you can build some other useful combinators.

--- goofy inefficient definitions
dup :: (Eq a, Eq b, BEnum a, BEnum b) => Rel a (a,a)
dup = rfan rid rid
swap ::(Eq a, Eq b, BEnum (a,b)) => Rel (a,b) (b,a)
swap = rfan rsnd rfst
par :: (Eq a, Eq c, BEnum a, BEnum c) => Rel a b -> Rel c d -> Rel (a,c) (b,d)
par f g =  rfan (rcompose f rfst) (rcompose g rsnd)

#### An Aside: Relations, Linear Algebra, Databases

The composition operation described above is not so unfamiliar as it may first appear.

Relation algebra has a great similarity to linear algebra. This connection can be made more clear by considering sparsity patterns of matrices and tensors. Sparsity patterns are a useful abstraction of linear algebraic operations. Instead of considering matrices of numbers, instead the entries are "zero" and "possibly nonzero" or, if you prefer, a matrix of boolean values corresponding to those same questions.

The ordinary row times column matrix multiplication corresponds to relation composition. Replace * with AND and + with OR. If any of the numbers is zero, then multiplying them will result in zero. In summing two numbers, if either is possibly nonzero, then the result is possibly nonzero.

Another interesting way of looking at it is that we are replacing the summation binding form $\sum_i$ with the logical quantifier $\exists_i$. Both introduce a scoped "dummy variable" i and have a lot of syntactic similarity. Other related forms include $\lambda i$, $\forall i$, $\int di$, $\max_i$ .

There is also an analog of the point free relation algebra in linear algebra. Linear algebra has the most widely used point free notation in the world, matrix notation. Consider the expressions $Ax=b$ and $X = ABC$ as compared to $\sum_j A_{ij} x_j = b_i$ and $X_{il} = \sum_{jk} A_{ij} B_{jk} C_{kl}$. Matrix notation is SO much better for certain calculations. Other pieces of the matrix notation include transpose, inverse, Kronecker product, the Khatri-Rao product, and Hadamard product. Their properties are more clear in the index free form in my opinion. I believe even massive tensor expressions can be written index free using these operators. There are also analogies to be drawn between the graphical notations in these different fields.

Databases can be thought of very similarly to sparse matrices. In principle, you could enumerate all the possible values for a column of a database. So you could think of a database as a giant matrix with a 1 if the item is in the database and 0 if not. Databases are very very sparse from this perspective, and you would never store them this way. The join operation is a relative of relational composition, however join usually operates via looking at the column names, whereas our join is position based.

Query optimization in databases has interesting analogs in sparse linear algebra. For example, the Taco compiler http://tensor-compiler.org/ is doing something very akin to a query optimizer.

#### Inverting Relations

Unlike functions, Relations are always "invertible". We call this the converse of a relation. When a function is invertible, it corresponds to the converse. In terms of the tuples underlying our representation, it just swaps them. Relations also possess operations trans and untrans that may be thought of as a kind of currying or as a partial inverse on a single parameter.

converse :: Rel a b -> Rel b a
converse = [(b,a) | (a,b) <- r]

untrans :: Rel a (b,c) -> Rel (a,b) c
untrans r = [((a,b),c)  | (a, (b,c)) <- r]

trans :: Rel (a,b) c -> Rel a (b, c)
trans r = [(a, (b,c))| ((a,b),c)  <- r]

Orderings can also be lifted to relations $(\leq) = \{(a,b) | a \leq b\}$. The composition of relations also respects the usual composition of ordering.

reflectOrd :: (Ord a, BEnum a) => Rel a a
reflectOrd = [(x,y) | x <- enumAll, y <- enumAll, x <= y]

Nondeterministic choice is sometimes represented in Haskell using Set returning functions a -> [b]. You may recall this from the context of the List monad. In fact in this case, we have an isomorphism as evidenced by tabulateSearch and searchRel.

tabulateSearch :: BEnum a => (a -> [b]) -> Rel a b
tabulateSearch f = [(a,b) | a <- enumAll, b <- f a]

searchRel :: Eq a => Rel a b -> (a -> [b])
searchRel r a = [b | (a', b) <- r, a == a']

Similarly partial functions can be reflected into relations

tabulatePartial :: BEnum a => (a -> Maybe b) -> Rel a b
tabulatePartial f = [(a,b) | a <- enumAll, b <- toList (f a)]

A useful trick is to lift sets/subsets to relations as a diagonal relation. $\{(a,a) | a \in S \}$. Projection onto the set can be achieve by composing with this relation. The identity results if you are talking about the entire set S.

diagRel :: [a] -> Rel a a
diagRel = map dup where dup x = (x,x)

leftSet :: Eq a => Rel a b -> [a]
leftSet = nub . (map fst)

rightSet :: Eq b => Rel a b -> [b]
rightSet = nub . (map snd)

#### Comparing Relations

We can compare sets by asking if one is a subset of the other $A\subseteq B$ . Relations can also be compared by this operation, which we call relational inclusion.

rSub :: (Eq a, Eq b) => Rel a b -> Rel a b -> Bool
rSub xs ys = and [x elem ys | x <- xs]

x <~ y = rSub x y

A subservient notion to this is relational equality.


rEq :: (Eq a, Eq b) => Rel a b -> Rel a b -> Bool
rEq xs ys = (xs rSub ys) && (ys rSub xs)

x ~~ y = rEq x y

Relational algebra is chockful of inequality style reasoning, which is richer and slightly more complicated than equality style reasoning. This is one of the benefits of moving from functional descriptions to a relational description.

Relations also form a lattice with respect to these comparisons. What the hell are lattices? In the context of finite relations, lattices may be over powered mathematical machinery, but it really is useful down the line. They give you binary operators that play nice with some kind of ordering, in our case relational inclusion. These two operations are the meet and the join, which find the greatest lower bound and least upper bound of the operands respectively. For our relations, these correspond to the more familiar notion of set intersection and union. The intersection of two sets is the biggest set that is in both of them. The union is the smallest set for which both sets are a subset of it.

meet' :: (Eq a, Eq b) => Rel a b -> Rel a b -> Rel a b
meet' xs ys = [x | x <- xs, x elem ys] -- intersection

join' :: (Eq a, Eq b) => Rel a b -> Rel a b -> Rel a b
join' p q = nub (p ++ q) -- union

instance (Eq a, Eq b) => Lattice (Rel a b) where
x /\ y = meet' x y
x \/ y = join' x y


Using meet/join vs intersection/union becomes more interesting when the domain is fancier than relations over finite domains. Issues of infinity can make this interesting, or when using a representation that can't explicitly represent arbitrary unions or intersections, but that instead has to approximate them. My favorite example is polyhedra. Polyhedra are not closed under unions. So in this case the join and union do not coincide. You need to take a convex hull of the union instead, which is the best approximation. Concretely, polyhedra can be represented as a list of their vertices, which generate the polyhedron. There is no way to express a union in this representation. Concatenating the lists represents taking the convex hull of the union.

An additional property that a lattice may possess is a largest and small element, called top ($\top$ ) and bottom ($\bot$). Our finite domain relations do have these.

instance (Eq a, Eq b, BEnum a, BEnum b) => BoundedMeetSemiLattice (Rel a b) where
top :: (BEnum a, BEnum b) => Rel a b
top = [(x,y) | x <- enumAll, y <- enumAll]
-- all possible tuples

-- because of superclass constraints :(
instance (Eq a, Eq b) => BoundedJoinSemiLattice (Rel a b) where
bottom :: Rel a b -- no tuples
bottom = [] 

#### Relational Division

And now finally we get to one of the most interesting, powerful, and confusing operations: relational division. Relational division is a kind of pseudo inverse to relational composition. In linear algebra, the pseudo inverse is a matrix that does the best job it can to invert another matrix in a least squares sense. If the matrix is actually invertible, it equals the inverse. Relational division does the best job it can to invert a relational composition. Instead of taking least squares as a criteria, it ensures that the result doesn't over approximate. If you had the inequality $X \cdot Y \subseteq Z$ and you want to solve for X, relational division is the thing that does that. The right division $Q = Z/Y$ is the largest relation such that $Q \cdot Y \subseteq Z$.

A helpful example is the similar operation of division in database tables.

And here is an implementation that I think is correct. I've goofed it up a couple times, it is a rather confusing construct.

rdiv :: (BEnum a, BEnum b, Eq a, Eq b, Eq c) => Rel a c -> Rel b c -> Rel a b
rdiv x y = [ (a,b)  | a <- enumAll, b <- enumAll, all (\c -> ((b,c) elem y)implies ((a,c) elem x)) (rightSet y)]

There also exists a very similar operation of ldiv.

Relational division encapsulates many notions of searching or optimizing. I invite you to read more about it in J.N. Oliveira's text or the Bird & de Moor text.

### Properties and QuickCheck

Relation algebra is so chock full of properties. This is a perfect opportunity for some QuickCheck , a randomized property testing framework. There are so many more to test. I need to dig through to collect up all the identities.

type R1 = Rel Bool Ordering

prop_ridleft :: Rel Bool Ordering -> Bool
prop_ridleft x = (rid <<< x) ~~ x

prop_ridright :: Rel Bool Ordering  -> Bool
prop_ridright x = (x <<< rid) ~~ x

prop_meet :: R1 -> R1  -> Bool
prop_meet x y = (x /\ y) <~ x

prop_meet' :: R1 -> R1  -> Bool
prop_meet' x y = (x /\ y) <~ y

prop_join_univ :: R1 -> R1 -> R1 -> Bool
prop_join_univ x y z = ((x \/ y) <~ z) == ((x <~ z) && (y <~ z))

prop_join :: R1 -> R1  -> Bool
prop_join x y = y <~ (x \/ y)

prop_meet_univ :: R1 -> R1 -> R1 -> Bool
prop_meet_univ x y z = (z <~ (x /\ y)) == ((z <~ x) && (z <~ y))

prop_top :: R1 -> Bool
prop_top x = x <~ top

prop_bottom :: R1 -> Bool
prop_bottom x = bottom <~ x

prop_bottom' :: R1 -> Bool
prop_bottom' x = (x <<< bottom) ~~ (bottom :: R1)

prop_trans_iso :: Rel (Bool, Ordering) Word8 -> Bool
prop_trans_iso x = untrans (trans x) == x

prop_rdiv :: Rel Bool Ordering -> Rel Word8 Ordering -> Bool
prop_rdiv g j = (j <<< (rdiv g j)) <~ g

prop_con :: R1 -> Bool
prop_con x = con (con x) ~~ x

prop_rdiv' :: Rel Bool Word8 -> Rel Bool Ordering -> Rel Word8 Ordering -> Bool
prop_rdiv' x g j = (x <~ (rdiv g j)) == ((j <<< x) <~ g)


### Bits and Bobbles

• Relations over continuous spaces. Vector subspaces (Linear Relations), Polyhedra (Linear inequality relations).
• Non Bool-valued Relations. Replace $\exists_x$ with $\max_x$. The weighted edgelist of a graph is a natural relation. By using composition we can ask about paths. We still have a comparison operator $\subseteq$ which now respects the ordering of weights
• Galois connections are cool.
• Relations combined with recursion schemes. Recursion schemes are the point free way of describing recursion.
• Moving into infinite spaces. How do we cope?
• Faster search. Some relations are best specified by functions, Maps, others, mixes and matching.
• If you "singletonize" relations a la the Agda project https://github.com/scmu/aopa, you get very interesting interconnections with profunctors, which people say are a categorical generalization of relations.
• Point-free DSLs are interesting and pleasant. Many worries about alpha renaming are gone, at the expense of point-free pain. A DSL like this may be necessary to choose good relational query plans

Edit: A follow up post on that type level angle here http://www.philipzucker.com/relational-algebra-with-fancy-types/

### References

Edit : A math exchange question about a -> [b] relational type. https://math.stackexchange.com/questions/3360026/can-division-be-expressed-intensionally-in-relation-algebra/3361351#3361351

Edit: An interesting comment and related library from /u/stevana

## Lens as a Divisibility Relation: Goofin’ Off With the Algebra of Types

Types have an algebra very analogous to the algebra of ordinary numbers (video). This is the basic table of correspondences. Code with all the extensions available here.

type a * b = (a,b)
type a + b = Either a b
type b ^ a = a -> b
type O = Void
type I = ()

-- check out typelits implementation

infixl 6 + -- ((a + b) + c) associates to the left
infixl 7 *
infixr 8 ^

-- derived definitions
type Succ a = I + a
type Two = Succ I
type Three = Succ Two
type Four = Succ Three

One way to see that this makes sense is by counting the cardinality of types built out of these combinators. Unit is the type with 1 inhabitant. Void has 0 inhabitants. If a has $n$ and b has $m$ possible values, then Either a b has $n + m$ inhabitants, (a,b) has $n*m$ and there are $n^m$ possible tabulations of a->b. We’re gonna stick to just polynomials for the rest of this, ignoring a->b.

Another way of looking at this is if two finitely inhabited types have the same number of inhabitants, then the types can be put into an isomorphism with each other. In other words, types modulo isomorphisms can be thought as representing the natural numbers. Because of this, we can build a curious proof system for the natural numbers using ordinary type manipulation.

In addition, we also get a natural way of expressing and manipulating polynomials.Polymorphic types can be seen as being very similar to polynomial expressions with natural coefficients N[x]. The polymorphic type variables have the ability to be instantiated to any value, corresponding to evaluating the polynomial for some number.

type ExamplePoly a = I + a + a * a * Three 

The Lens ecosystem gives some interesting combinators for manipulating this algebra. The type Iso' a b contains isomorphisms. Since we’re only considering types up to isomorphism, this Iso' represents equality. We can give identity isomorphisms, compose isomorphisms and reverse isomorphisms.

type a ~~ b = Iso' a b

refl ::  a ~~ a
refl = id

compose :: (a ~~ b) -> (b ~~ c) -> (a ~~ c)
compose = (.)

rev :: a ~~ b -> b ~~ a
rev = from

We can already form a very simple proof.

-- a very simple proof. Holds by definition
oneonetwo :: (I + I) ~~ Two
oneonetwo = id

Now we’ll add some more combinators, basically the axioms that the types mod isos are a commutative semiring. Semirings have an addition and multiplication operator that distribute over each other. It is interesting to note that I believe all of these Iso' actually are guaranteed to be isomorphisms ( to . from = id and from . to = id ) because of parametricity. from and to are unique ignoring any issues with bottoms because the polymorphic type signature is so constraining. This is not usually guaranteed to be true in Haskell just from saying it is an Iso'. If I give you an Iso' Bool Bool it might actually be the iso (const True) (const True) for example, which is not an isomorphism.

-- now we start to state our axioms
-- interestingly, I believe the Iso and Lens laws to follow actually follow via parametricity.

-- associativity + a identity object make mul and plus a monoid
plus_assoc :: Iso' (a + (b + c)) ((a + b) + c)
plus_assoc =  iso assoc unassoc  where
assoc (Left a) = Left (Left a)
assoc (Right (Left b)) = Left (Right b)
assoc (Right (Right c)) = Right c
unassoc (Left (Left a)) = Left a
unassoc (Left (Right b)) = Right (Left b)
unassoc (Right c) = (Right (Right c))

mul_assoc :: Iso' (a * (b * c)) ((a * b) * c)
mul_assoc =  iso (\(a,(b,c)) -> ((a,b),c)) (\((a,b),c) -> (a,(b,c)))

-- use absurd from Data.Void for empty case.
id_plus :: Iso' (a + O) a
id_plus = iso (\case Left x -> x
Right x -> absurd x) Left

id_mul :: Iso' (a * I) a
id_mul = iso (\(x,_) -> x) (\x -> (x,()))

-- they are also commutative
-- specialized version of swapped from Control.Lens.Iso for future clarity
comm_plus :: Iso' (a + b) (b + a)
comm_plus = swapped
comm_mul :: Iso' (a * b) (b * a)
comm_mul = swapped

-- I don't see this one in Lens. Maybe it is there though?
-- distributive property
rdist :: Iso' (a * (b + c)) (a * b + a * c)
rdist = iso (\(a,bc) -> (a,) +++ (a,) \$ bc) (\case Left (a,b) -> (a, Left b)
Right (a,c) -> (a, Right c))

mul_zero :: Iso' (a,O) O
mul_zero = iso (\(_,y) -> absurd y) absurd


There are also combinators for lifting isomorphisms into bifunctors: firsting, seconding, and bimapping. These are important for indexing into subexpressions of our types in a point-free style.

-- a specialized version of firsting and seconding for clarity
-- firsting and seconding feel to me more like words for tuples
lefting :: (a ~~ b) -> (a + c) ~~ (b + c)
lefting = firsting

righting :: (a ~~ b) -> (c + a) ~~ (c + b)
righting = seconding

Here is a slightly more complicated proof now available to us.

-- a more complicated proof
twotwofour :: Iso' (Two + Two) Four
twotwofour = rev plus_assoc

We can attempt a more interesting and difficult proof. I developed this iteratively using . _ hole expressions so that GHC would tell me what I had manipulated my type to at that point in my proof.

ldist ::   ((b + c) * a) ~~ (b * a + c * a)
ldist = comm_mul . rdist . (lefting comm_mul) . (righting comm_mul)

-- very painful. Using holes _ and error messages absolutely essential
factorexample ::  ((a + I) * (a + I)) ~~ (a * a + Two * a + I)
factorexample = rdist .  -- distribute out the right side
(lefting (comm_mul . rdist)) . -- in the left sum term distribute out
(righting (comm_mul . rdist)) . -- in the right sum term distribute out
plus_assoc . -- reassociate plus to the left
(lefting (lefting (righting comm_mul))) . -- turn a * I term to I * a
(lefting (rev plus_assoc)) . -- associate the two a*I terms together into an (a * I + a * I) term
(lefting (righting (rev ldist))) . -- factor together that subterm into Two * a
(righting id_mul) -- convert last I * I term into I



The proof here is actually pretty trivial and can be completely automated away. We’ll get to that later.

If Iso' is equality, what about the other members of the Lens family? Justin Le says that Lens s a are witness to the isomorphism of a type s to the tuple of something and a. Prism witness a similar thing for sums. Once we are only considering types mod isos, if you think about it, these are expressions of two familiar relations on the natural numbers: the inequality relation and the divisibility relation

type a >= b = Prism' a b
type a .| b = Lens' a b

Mathematically, these relations can be composed with equalities, just like in the lens hierarchy Lens and Prism can be composed with Iso. Both form a category, since they both have id and (.).

{-

the core combinators from the lens library for manipulating these are

_1 :: (a * b) .| a
_2 :: (a * b) .| b

_Left :: (a + b) >= a
_Right :: (a + b) >= b

For example:
-}
twodiv4 :: (Two * Two) .| Two
twodiv4 = _1

onelesstwo :: Two >= I
onelesstwo = _Left

threedivthree :: Three .| Three
threedivthree = id


Here are a couple identities that we can’t derive from these basic combinators. There are probably others. Woah-ho, my bad. These are totally derivable using id_mul, id_plus, mul_zero, _1, _2, _Left, _Right.

-- once again, these are true via parametricity, more or less
one_div :: a .| I
-- one_div = lens (const ()) const
one_div = (rev id_mul) . _2

zero_lte :: a >= O
-- zero_lte = prism' absurd (const Nothing)
zero_lte = (rev id_plus) . _Right

zero_div :: O .| a
-- zero_div = lens absurd const
zero_div = (rev mul_zero) . _1

Pretty neat! Random thoughts and questions before we get into the slog of automation:

• Traversal is the “is polynomial in” relation, which seems a rather weak statement on it’s own.
• Implementing automatic polynomial division is totally possible and interesting
• What is the deal with infinite types like [a]? Fix. I suppose this is a theory of formal power series / combinatorial species. Fun combinatorics, generatingfunctionology. Brent Yorgey did his dissertation on this stuff. Wow. I’ve never really read this. It is way more relevant than I realized.
• Multivariate polynomial algorithms would also be interesting to explore (Grobner basis, multivariate division)
• Derivatives of types and zippers – Conor McBride
• Negative Numbers and Fractions?
• Lifting to rank-1 types. Define Negative and Fractions via Galois connection?
-- partially applied +
type P n = (Either n)
-- partially applied *
type M n = (,) n

Edit: /u/carette (wonder who that is ðŸ˜‰ ) writes:

“You should dig into
J Carette, A Sabry Computing with semirings and weak rig groupoids, in Proceedings of ESOP 2016, p. 123-148. Agda code in https://github.com/JacquesCarette/pi-dual/tree/master/Univalence. A lot of the algebra you develop is there too.

If you hunt around in my repos, you’ll also find things about lenses, exploring some of the same things you mention here.”

Similar ideas taken further and with more sophistication. Very interesting. Check it out.

## Automation

Our factor example above was quite painful, yet the theorem was exceedingly obvious by expansion of the left and right sides. Can we automate that? Yes we can!

Here’s the battle plan:

• Distribute out all expressions like $a*(b+c)$ so that all multiplication nodes appear at the bottom of the tree.
• Reduce the expression by absorbing all stupid $a*1$, $a*0$, $a+0$ terms.
• Reassociate everything to the right, giving a list like format
• Sort the multiplicative terms by power of the variable

Once we have these operations, we’ll combine them into a canonical operation. From there, most equality proofs can be performed using the rewrite operation, which merely puts both sides into canonical form

-- | The main combinator to expand out any polynomial into canonical form with terms sorted

canonical :: (Dist a b, Absorb b c, RightAssoc c d, SortTerm d e) => a ~~ e
canonical = dist . absorb . rightAssoc . sortTerm

rewrite :: forall a' a b c d e e' d' c' b'. (Dist a b, Absorb b c, RightAssoc c d,
SortTerm d e, e ~ e', Dist a' b',
Absorb b' c', RightAssoc c' d', SortTerm d' e') => a ~~ a'
rewrite = canonical . (rev canonical)

Once we have those, the painful theorem above and harder ones becomes trivial.

ex9 :: (a ~ V a') => ((a + I) * (I + a)) ~~ (a * a + Two * a + I)
ex9 = rewrite

ex10 :: (a ~ V a') => ((a + I) * (I + a) * (Two + a)) ~~  (a * a * a + a * a * Four + Two + Four * a + a)
ex10 = rewrite 

Now we’ll build the Typeclasses necessary to achieve each of these aims in order. The Typeclass system is perfect for what we want to do, as it builds terms by inspecting types. It isn’t perfect in the sense that typeclass pattern matching needs to be tricked into doing what we need. I have traded in cleverness and elegance with verbosity.

In order to make our lives easier, we’ll need to tag every variable name with a newtype wrapper. Otherwise we won’t know when we’ve hit a leaf node that is a variable. I’ve used this trick before here in an early version of my faking Compiling to Categories series. These wrappers are easily automatically stripped.

-- For working with Variables

-- a newtype to mark variable position
newtype V a = V a

-- Suggested usage, bind the V in a unification predicate to keep expression looking clean
-- (V a' ~ a) => (a + I) * a

-- multinomials. to be implemented some other day
-- a phantom l labelled newtype for variable ordering.
newtype VL l a = VL a

A common pattern I exploit is to use a type family to drive complicated recursion. Closed type families allow more overlap and default patterns which is very useful for programming. However, type families do not carry values, so we need to flip flop between the typeclass and type family system to achieve our ends.

Here is the implementation of the distributor Dist. We make RDist and LDist typeclasses that make a sweep of the entire tree, using ldist and rdist as makes sense. It was convenient to separate these into two classes for my mental sanity. I am not convinced even now that I have every case. Then the master control class Dist runs these classes until any node that has a (*) in it has no nodes with (+) underneath, as checked by the HasPlus type family.


class RDist a b | a -> b where
rdist' :: Iso' a b

instance RDist a a' => RDist (a * I) (a' * I) where
rdist' = firsting rdist'
instance RDist a a' => RDist (a * O) (a' * O) where
rdist' = firsting rdist'
instance RDist a a' => RDist (a * (V b)) (a' * (V b)) where
rdist' = firsting rdist'

instance (RDist (a * b) ab, RDist (a * c) ac) => RDist (a * (b + c)) (ab + ac) where
rdist' = rdist . (bimapping rdist' rdist')
instance (RDist a a', RDist (b * c) bc) => RDist (a * (b * c)) (a' * bc) where
rdist' = (bimapping rdist' rdist')
instance (RDist a a', RDist b b') => RDist (a + b) (a' + b') where
rdist' = bimapping rdist' rdist'

instance RDist O O where
rdist' = id
instance RDist I I where
rdist' = id
instance RDist (V a) (V a) where
rdist' = id

-- can derive ldist from swapped . rdist' . swapped?

class LDist a b | a -> b where
ldist' :: Iso' a b
instance (LDist (b * a) ab, LDist (c * a) ac) => LDist ((b + c) * a) (ab + ac) where
ldist' = ldist . (bimapping ldist' ldist')
instance (LDist (b * c) bc, LDist a a') => LDist ((b * c) * a) (bc * a') where
ldist' = bimapping ldist' ldist'
instance LDist a a' => LDist (I * a) (I * a') where
ldist' = seconding ldist'
instance LDist a a' => LDist (O * a) (O * a') where
ldist' = seconding ldist'
instance LDist a a' => LDist ((V b) * a) ((V b) * a') where
ldist' = seconding ldist'

instance (LDist a a', LDist b b') => LDist (a + b) (a' + b') where
ldist' = bimapping ldist' ldist'

instance LDist O O where
ldist' = id
instance LDist I I where
ldist' = id
instance LDist (V a) (V a) where
ldist' = id

type family HasPlus a where
HasPlus (a + b) = 'True
HasPlus (a * b) = (HasPlus a) || (HasPlus b)
HasPlus I = 'False
HasPlus O = 'False
HasPlus (V _) = 'False

class Dist' f a b | f a -> b where
dist' :: a ~~ b
instance (f ~ HasPlus ab', LDist (a * b) ab,
RDist ab ab', Dist' f ab' ab'') => Dist' 'True (a * b) ab'' where
dist' = ldist' . rdist' . (dist' @f)
instance Dist' 'False (a * b) (a * b) where
dist' = id
instance (HasPlus a ~ fa, HasPlus b ~ fb,
Dist' fa a a', Dist' fb b b') => Dist' x (a + b) (a' + b') where
dist' = bimapping (dist' @fa) (dist' @fb)
-- is that enough though? only dist if
instance Dist' x I I where
dist' = id
instance Dist' x O O where
dist' = id
instance Dist' x (V a) (V a) where
dist' = id

class Dist a b | a -> b where
dist :: a ~~ b
-- dist' should distributed out all multiplications
instance (f ~ HasPlus a, Dist' f a b) => Dist a b where
dist = dist' @f


Next is the Absorb type class. It is arranged somewhat similarly to the above. Greedily absorb, and keep doing it until no absorptions are left. I think that works.


-- RAbsorb matches only on the right hand side of binary operators
-- matching on both sides is ungainly to write
class RAbsorb a b | a -> b where
rabsorb :: a ~~ b
-- obvious absorptions
instance RAbsorb x x' => RAbsorb (x + O) x' where
rabsorb = id_plus . rabsorb
instance RAbsorb x x' => RAbsorb (x + I) (x' + I) where
rabsorb = firsting rabsorb
instance RAbsorb x x' => RAbsorb (x * I) x' where
rabsorb = id_mul . rabsorb
instance RAbsorb (x * O) O where
rabsorb = mul_zero
instance RAbsorb x x' => RAbsorb (x * (V a)) (x' * (V a)) where
rabsorb = firsting rabsorb
instance RAbsorb x x' => RAbsorb (x + (V a)) (x' + (V a)) where
rabsorb = lefting rabsorb
-- recursion steps
instance (RAbsorb (y * z) yz, RAbsorb x x') => RAbsorb (x * (y * z)) (x' * yz) where
rabsorb = bimapping rabsorb rabsorb
instance (RAbsorb (y + z) yz, RAbsorb x x') => RAbsorb (x * (y + z)) (x' * yz) where
rabsorb = bimapping rabsorb rabsorb
instance (RAbsorb (y + z) yz, RAbsorb x x') => RAbsorb (x + (y + z)) (x' + yz) where
rabsorb = bimapping rabsorb rabsorb
instance (RAbsorb (y * z) yz, RAbsorb x x') => RAbsorb (x + (y * z)) (x' + yz) where
rabsorb = bimapping rabsorb rabsorb
-- base cases
instance RAbsorb O O where
rabsorb = id
instance RAbsorb I I where
rabsorb = id
instance RAbsorb (V a) (V a) where
rabsorb = id

-- mirror of RAbsorb
class LAbsorb a b | a -> b where
labsorb :: a ~~ b
instance LAbsorb x x' => LAbsorb (O + x) x' where
labsorb = comm_plus . id_plus . labsorb
instance LAbsorb x x' => LAbsorb (I + x) (I + x') where
labsorb = righting labsorb
instance LAbsorb x x' => LAbsorb (I * x) x' where
labsorb = comm_mul . id_mul . labsorb
instance LAbsorb (O * x) O where
labsorb = comm_mul . mul_zero
instance LAbsorb x x' => LAbsorb ((V a) + x) ((V a) + x') where
labsorb = righting labsorb
instance LAbsorb x x' => LAbsorb ((V a) * x) ((V a) * x') where
labsorb = seconding labsorb

instance (LAbsorb (y * z) yz, LAbsorb x x') => LAbsorb ((y * z) * x) (yz * x') where
labsorb = bimapping labsorb labsorb
instance (LAbsorb (y + z) yz, LAbsorb x x') => LAbsorb ((y + z) * x) (yz * x') where
labsorb = bimapping labsorb labsorb
instance (LAbsorb (y + z) yz, LAbsorb x x') => LAbsorb ((y + z) + x) (yz + x') where
labsorb = bimapping labsorb labsorb
instance (LAbsorb (y * z) yz, LAbsorb x x') => LAbsorb ((y * z) + x) (yz + x') where
labsorb = bimapping labsorb labsorb

instance LAbsorb O O where
labsorb = id
instance LAbsorb I I where
labsorb = id
instance LAbsorb (V a) (V a) where
labsorb = id

-- labsorb :: (Swapped p, RAbsorb (p b a) (p b' a')) => (p a b) ~~ (p a' b')
-- labsorb = swapped . rabsorb . swapped

-- a test function to see if an expression is properly absorbed
type family Absorbed a where
Absorbed (O + a) = 'False
Absorbed (a + O) = 'False
Absorbed (a * I) = 'False
Absorbed (I * a) = 'False
Absorbed (a * O) = 'False
Absorbed (O * a) = 'False
Absorbed (a + b) = (Absorbed a) && (Absorbed b)
Absorbed (a * b) = (Absorbed a) && (Absorbed b)
Absorbed I = 'True
Absorbed O = 'True
Absorbed (V a) = 'True

-- iteratively rabsorbs and leftabsorbs until tree is properly absorbed.
class Absorb' f a b | f a -> b where
absorb' :: a ~~ b
instance Absorb' 'True a a where
absorb' = id
instance (LAbsorb a a', RAbsorb a' a'',
f ~ Absorbed a'', Absorb' f a'' a''') => Absorb' 'False a a''' where
absorb' = labsorb . rabsorb . (absorb' @f)

-- wrapper class to avoid the flag.
class Absorb a b | a -> b where
absorb :: a ~~ b
instance (f ~ Absorbed a, Absorb' f a b) => Absorb a b where
absorb = absorb' @f


The Associators are a little simpler. You basically just look for the wrong association pattern and call plus_assoc or mul_assoc until they don’t occur anymore, then recurse. We can be assured we’re always making progress if we either switch some association structure or recurse into subparts.


class LeftAssoc a b | a -> b where
leftAssoc :: Iso' a b

instance LeftAssoc a a' => LeftAssoc (a + I) (a' + I) where
leftAssoc = firsting leftAssoc
instance LeftAssoc a a' => LeftAssoc (a + O) (a' + O) where
leftAssoc = firsting leftAssoc
instance LeftAssoc a a' => LeftAssoc (a * I) (a' * I) where
leftAssoc = firsting leftAssoc
instance LeftAssoc a a' => LeftAssoc (a * O) (a' * O) where
leftAssoc = firsting leftAssoc
instance LeftAssoc a a' => LeftAssoc (a * (V b)) (a' * (V b)) where
leftAssoc = firsting leftAssoc
instance LeftAssoc a a' => LeftAssoc (a + (V b)) (a' + (V b)) where
leftAssoc = firsting leftAssoc

instance (LeftAssoc ((a + b) + c) abc') => LeftAssoc (a + (b + c)) abc' where
leftAssoc = plus_assoc . leftAssoc
instance (LeftAssoc ((a * b) * c) abc') => LeftAssoc (a * (b * c)) abc' where
leftAssoc = mul_assoc . leftAssoc

instance (LeftAssoc (b * c) bc, LeftAssoc a a') => LeftAssoc (a + (b * c)) (a' + bc) where
leftAssoc = bimapping leftAssoc leftAssoc
-- a * (b + c) ->  a * b + a * c
-- This case won't happen if we've already distribute out.
instance (LeftAssoc (b + c) bc, LeftAssoc a a') => LeftAssoc (a * (b + c)) (a' * bc) where
leftAssoc = bimapping leftAssoc leftAssoc

instance LeftAssoc O O where
leftAssoc = id
instance LeftAssoc I I where
leftAssoc = id
instance LeftAssoc (V a) (V a) where
leftAssoc = id

-- right assoc will completely associate strings of + or -. Mixed terms are not associated.
-- cases  on left hand side of binary expression
-- always makes progress by either reassociating or recursing
class RightAssoc a b | a -> b where
rightAssoc :: Iso' a b
instance (RightAssoc (a + (b + c)) abc') => RightAssoc ((a + b) + c) abc' where
rightAssoc = (rev plus_assoc) . rightAssoc
instance (RightAssoc (a * (b * c)) abc') => RightAssoc ((a * b) * c) abc' where
rightAssoc = (rev mul_assoc) . rightAssoc
instance RightAssoc a a' => RightAssoc (I + a) (I + a') where
rightAssoc = seconding rightAssoc
instance RightAssoc a a' => RightAssoc (O + a) (O + a') where
rightAssoc = seconding rightAssoc
instance RightAssoc a a' => RightAssoc (I * a) (I * a') where
rightAssoc = seconding rightAssoc
instance RightAssoc a a' => RightAssoc (O * a) (O * a') where
rightAssoc = seconding rightAssoc
instance RightAssoc a a' => RightAssoc ((V b) + a) ((V b) + a') where
rightAssoc = seconding rightAssoc
instance RightAssoc a a' => RightAssoc ((V b) * a) ((V b) * a') where
rightAssoc = seconding rightAssoc

instance (RightAssoc (b * c) bc, RightAssoc a a') => RightAssoc ((b * c) + a) (bc + a') where
rightAssoc = bimapping rightAssoc rightAssoc
instance (RightAssoc (b + c) bc, RightAssoc a a') => RightAssoc ((b + c) * a) (bc * a') where
rightAssoc = bimapping rightAssoc rightAssoc

instance RightAssoc O O where
rightAssoc = id
instance RightAssoc I I where
rightAssoc = id
instance RightAssoc (V a) (V a) where
rightAssoc = id


Finally, the SortTerm routine. SortTerm is a bubble sort. The typeclass Bubble does a single sweep of swapping down the type level list-like structure we’ve built. The SortTerm uses the Sorted type family to check if it is finished. If it isn’t, it call Bubble again.


type family (a :: k) == (b :: k) :: Bool where
f a == g b = f == g && a == b
a == a = 'True
_ == _ = 'False

type family SortedTerm a :: Bool where
SortedTerm (a + (b + c)) = (((CmpTerm a b) == 'EQ) || ((CmpTerm a b) == 'GT)) && (SortedTerm (b + c))
SortedTerm (a + b) = ((CmpTerm a b) == 'EQ) || ((CmpTerm a b) == 'GT)
--SortedTerm a = 'True
SortedTerm I = 'True
SortedTerm O = 'True
SortedTerm (V a) = 'True

-- higher powers of V are bigger.
-- CmpTerm compares TimesLists.
type family CmpTerm a b where
CmpTerm ((V a) * b) ((V a) * c) = CmpTerm b c
CmpTerm ((V a) * b) (V a) = 'GT
CmpTerm (V a) ((V a) * b)  = 'LT
CmpTerm I ((V a) * b) = 'LT
CmpTerm ((V a) * b) I = 'GT
CmpTerm (V a) (V a) = 'EQ
CmpTerm I (V a) = 'LT
CmpTerm (V a) I = 'GT
CmpTerm I I = 'EQ

-- Maybe this is all uneccessary since we'll expand out and abosrb to a*a + a*a + a + a + a + a

-- type a == b = TEq.(==) a b

-- Head and Tail of PlusLists
PlusHead (x + y) = x
type family PlusTail a where
PlusTail (x + y) = y

-- bubble assume a plusList of multiplicative terms. I.E. fully distributed, fully rightassociated , fully absorbed
-- does one pass of a bubble sort
class Bubble f a b | f a -> b where
bubble :: a ~~ b
-- more to go
instance (f ~ CmpTerm b (PlusHead c), Bubble f (b + c) bc) => Bubble 'EQ (a + (b + c)) (a + bc) where
bubble = righting (bubble @f)
instance (f ~ CmpTerm b (PlusHead c), Bubble f (b + c) bc) => Bubble 'GT (a + (b + c)) (a + bc) where
bubble = righting (bubble @f)
instance (f ~ CmpTerm a (PlusHead c), Bubble f (a + c) ac) => Bubble 'LT (a + (b + c)) (b + ac) where
bubble = plus_assoc . (lefting comm_plus) . (rev plus_assoc) . righting (bubble @f)

-- The times, or constants shows that we're at the end of our + list.
instance Bubble 'EQ (a + (b * c)) (a + (b * c)) where
bubble = id
instance Bubble 'GT (a + (b * c)) (a + (b * c)) where
bubble = id
instance Bubble 'LT (a + (b * c)) ((b * c) + a) where
bubble = comm_plus

instance Bubble 'EQ (a + I) (a + I) where
bubble = id
instance Bubble 'GT (a + I) (a + I) where
bubble = id
instance Bubble 'LT (a + I) (I + a) where
bubble = comm_plus

instance Bubble 'EQ (a + O) (a + O) where
bubble = id
instance Bubble 'GT (a + O) (a + O) where
bubble = id
instance Bubble 'LT (a + O) (O + a) where -- shouldn't happen
bubble = comm_plus

instance Bubble 'EQ (a + (V b)) (a + (V b)) where
bubble = id
instance Bubble 'GT (a + (V b)) (a + (V b)) where
bubble = id
instance Bubble 'LT (a + (V b)) ((V b) + a) where
bubble = comm_plus
-- goofy base cases in case bubble gets called on a single element
instance Bubble x O O where
bubble = id
instance Bubble x I I where
bubble = id
instance Bubble x (V a) (V a) where
bubble = id

-- sort term assumes rightassociated, fully distributed, fully I O absorbed expressions
class SortTerm' f a b | f a -> b where -- f is flag whether PlusList is sorted.
sortTerm' :: a ~~ b
instance SortTerm' 'True a a where
sortTerm' = id
-- a single term with no plus shouldn't get here. That is why PlusTail is ok.
instance (f ~ CmpTerm (PlusHead a) (PlusHead (PlusTail a)), Bubble f a a',
f' ~ SortedTerm a', SortTerm' f' a' b) => SortTerm' 'False a b where
sortTerm' = (bubble @f) . (sortTerm' @f')
class SortTerm a b | a -> b where -- f is flag whether PlusList is sorted.
sortTerm :: a ~~ b
instance (f ~ SortedTerm a, SortTerm' f a a') => SortTerm a a' where
sortTerm = sortTerm' @f


Hope you thought this was neat!

## Annihilation Creation with Wick Contraction in Python

Trying an alternative approach to dealing with the algebra of qft computationally based on wick contraction. We can easily find an algorithm that finds all possible pairs. Then we just reduce it all accordingly. Some problems: The combinatorics are going to freak out pretty fast.

I think my use of functional stuff like maps and lambdas is intensely unclarifying the code.

import numpy as np

#Needs to return a list of lists of pairs
def pair(mylist):
#Base case
if len(mylist) == 2:
return [[(mylist[0], mylist[1])]]
#popoff the first element
element1 = mylist[0]
pairs = []
for i in range(1,len(mylist)):
#Pick one
element2 = mylist[i]
#get all the pairs expluidng the already picked pair
subpairs = pair(mylist[1:i] + mylist[i+1:])
#Put the picked pair back in x is a list of pairs.
pairs = pairs + map(lambda x: [(element1,element2)] + x ,subpairs)
return pairs

#For Fermionic pairing, it might be conveneitn to extend pairs to (element1,element2,+-1)
#With the sign depending on wheter i is even or odd

'''
pairing = pair([1,2,3,4])
print pairing
print len(pairing)
'''

#Here's an idea, I can use a and adag as objects with indices
#and return functions (two point green's functions)
# or could return matrcies that are discretized green's functions.
def contract(twoguys):
if twoguys[0] == 'a' and twoguys[1] == 'adag':
return 1
else:
return 0

#Confusing, but what