Linear Relation Algebra of Circuits with HMatrix

Oooh this is a fun one.

I’ve talked before about relation algebra and I think it is pretty neat. http://www.philipzucker.com/a-short-skinny-on-relations-towards-the-algebra-of-programming/. In that blog post, I used finite relations. In principle, they are simple to work with. We can perform relation algebra operations like composition, meet, and join by brute force enumeration.

Unfortunately, brute force may not always be an option. First off, the finite relations grow so enormous as to be make this infeasible. Secondly, it is not insane to talk about relations or regions with an infinite number of elements, such as some continuous blob in 2D space. In that case, we can’t even in principle enumerate all the points in the region. What are we to do? We need to develop some kind of finite parametrization of regions to manipulate. This parametrization basically can’t possibly be complete in some sense, and we may choose more or less powerful systems of description for computational reasons.

In this post, we are going to be talking about linear or affine subspaces of a continuous space. These subspaces are hyperplanes. Linear subspaces have to go through the origin, while affine spaces can have an offset from the origin.

In the previous post, I mentioned that the finite relations formed a lattice, with operations meet and join. These operations were the same as set intersection and union so the introduction of the extra terminology meet and join felt a bit unwarranted. Now the meet and join aren’t union and intersection anymore. We have chosen to not have the capability to represent the union of two vectors, instead we can only represent the smallest subspace that contains them both, which is the union closed under vector addition. For example, the join of a line and point will be the plane that goes through both.

Linear/Affine stuff is great because it is so computational. Most questions you cant to ask are answerable by readily available numerical linear algebra packages. In this case, we’ll use the Haskell package HMatrix, which is something like a numpy/scipy equivalent for Haskell. We’re going to use type-level indices to denote the sizes and partitioning of these spaces so we’ll need some helper functions.

Matrices are my patronum. They make everything good. Artwork courtesy of David

In case I miss any extensions, make typos, etc, you can find a complete compiling version here https://github.com/philzook58/ConvexCat/blob/master/src/LinRel.hs

type BEnum a = (Enum a, Bounded a) 

-- cardinality. `size` was already taken by HMatrix :(
card :: forall a. (BEnum a) => Int
card = (fromEnum (maxBound @a)) - (fromEnum (minBound @a)) + 1

In analogy with sets of tuples for defining finite relations, we partition the components of the linear spaces to be “input” and “output” indices/variables \begin{bmatrix} x_1 & x_2 & x_3 & ... & y_1 & y_2 & y_3 & ... \end{bmatrix}. This partition is somewhat arbitrary and easily moved around, but the weakening of strict notions of input and output as compared to functions is the source of the greater descriptive power of relations.

Relations are extensions of functions, so linear relations are an extension of linear maps. A linear map has the form y = Ax. A linear relation has the form Ax + By = 0. An affine map has the form y = Ax + b and an affine relation has the form Ax + By = b.

There are at least two useful concrete representation for subspaces.

  1. We can write a matrix A and vector b down that corresponds to affine constraints. Ax = b. The subspace described is the nullspace of A plus a solution of the equation. The rows of A are orthogonal to the space.
  2. We can hold onto generators of subspace. x = A' l+b where l parametrizes the subspace. In other words, the subspace is generated by / is the span of the columns of A'. It is the range of A'.

We’ll call these two representations the H-Rep and V-Rep, borrowing terminology from similar representations in polytopes (describing a polytope by the inequalities that define it’s faces or as the convex combination of it’s vertices). https://inf.ethz.ch/personal/fukudak/lect/pclect/notes2015/PolyComp2015.pdf These two representations are dual in many respects.

-- HLinRel holds A x = b constraint
data HLinRel a b = HLinRel (Matrix Double) (Vector Double) deriving Show

-- x = A l + b. Generator constraint. 
data VLinRel a b = VLinRel (Matrix Double) (Vector Double) deriving Show

It is useful to have both reps and interconversion routines, because different operations are easy in the two representations. Any operations defined on one can be defined on the other by sandwiching between these conversion functions. Hence, we basically only need to define operations for one of the reps (if we don’t care too much about efficiency loss which, fair warning, is out the window for today). The bulk of computation will actually be performed by these interconversion routines. The HMatrix function nullspace performs an SVD under the hood and gathers up the space with 0 singular values.

-- if A x = b then x is in the nullspace + a vector b' solves the equation
h2v :: HLinRel a b -> VLinRel a b
h2v (HLinRel a b) = VLinRel a' b' where
        b' = a <\> b -- least squares solution
        a' = nullspace a

-- if x = A l + b, then A' . x = A' A l + A' b = A' b because A' A = 0
v2h :: VLinRel a b -> HLinRel a b
v2h (VLinRel a' b') = HLinRel a b where
        b = a #> b' -- matrix multiply
        a = tr $ nullspace (tr a') -- orthogonal space to range of a.
-- tr is transpose and not trace? A little bit odd, HMatrix.

These linear relations form a category. I’m not using the Category typeclass because I need BEnum constraints hanging around. The identity relations is x = y aka Ix - Iy = 0.

hid :: forall a. BEnum a => HLinRel a a
hid =  HLinRel (i ||| (- i)) (vzero s) where 
                            s = card @a
                            i = ident s

Composing relations is done by combining the constraints of the two relations and then projecting out the interior variables. Taking the conjunction of constraints is easiest in the H-Rep, where we just need to vertically stack the individual constraints. Projection easily done in the V-rep, where you just need to drop the appropriate section of the generator vectors. So we implement this operation by flipping between the two.

hcompose :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel b c -> HLinRel a b -> HLinRel a c
hcompose (HLinRel m b) (HLinRel m' b') = let a'' = fromBlocks [[       ma',           mb' ,    0       ],
                                                               [         0 ,    mb,        mc          ]] in
                                         let b'' = vjoin [b', b] in 
                                         let (VLinRel q p) = h2v (HLinRel a'' b'') in -- kind of a misuse
                                         let q' = (takeRows ca q)  -- drop rows belonging to @b
                                                       === 
                                                  (dropRows (ca + cb) q) in
                                         let [x,y,z] =  takesV [ca,cb,cc] p in
                                         let p'=  vjoin [x,z] in -- rebuild without rows for @b
                                         v2h (VLinRel q' p') -- reconstruct HLinRel
                                       where 
                                           ca = card @a
                                           cb = card @b 
                                           cc = card @c
                                           sb = size b -- number of constraints in first relation
                                           sb' = size b' -- number of constraints in second relation
                                           ma' = takeColumns ca m'
                                           mb' = dropColumns ca m'
                                           mb = takeColumns cb m
                                           mc = dropColumns cb m

(<<<) :: forall a b c. (BEnum a, BEnum b, BEnum c) => HLinRel b c -> HLinRel a b -> HLinRel a c
(<<<) = hcompose

We can implement the general cadre of relation operators, meet, join, converse. I feel the converse is the most relational thing of all. It makes inverting a function nearly a no-op.

hjoin :: HLinRel a b -> HLinRel a b -> HLinRel a b
hjoin v w = v2h $ vjoin' (h2v v) (h2v w)

-- hmatrix took vjoin from me :(
-- joining means combining generators and adding a new generator
-- Closed under affine combination l * x1 + (1 - l) * x2 
vjoin' :: VLinRel a b -> VLinRel a b -> VLinRel a b
vjoin' (VLinRel a b) (VLinRel a' b') = VLinRel (a ||| a' ||| (asColumn (b - b'))) b

-- no constraints, everything
-- trivially true
htop :: forall a b. (BEnum a, BEnum b) => HLinRel a b 
htop = HLinRel (vzero (1,ca + cb)) (konst 0 1) where 
                                      ca = card @a
                                      cb = card @b 

-- hbottom?
                              
hconverse :: forall a b. (BEnum a, BEnum b) => HLinRel a b -> HLinRel b a 
hconverse (HLinRel a b) = HLinRel ( (dropColumns ca a) |||  (takeColumns ca a)) b where 
    ca = card @a
    cb = card @b  

Relational inclusion is the question of subspace inclusion. It is fairly easy to check if a VRep is in an HRep (just see plug the generators into the constraints and see if they obey them) and by using the conversion functions we can define it for arbitrary combos of H and V.


-- forall l. A' ( A l + b) == b'
-- is this numerically ok? I'm open to suggestions.
vhsub :: VLinRel a b -> HLinRel a b -> Bool
vhsub (VLinRel a b) (HLinRel a' b') = (naa' <=  1e-10 * (norm_2 a') * (norm_2 a)  ) && ((norm_2 ((a' #> b) - b')) <= 1e-10 * (norm_2 b')  ) where
          naa' = norm_2 (a' <> a)

hsub :: HLinRel a b -> HLinRel a b -> Bool
hsub h1 h2 = vhsub (h2v h1) h2

heq :: HLinRel a b -> HLinRel a b -> Bool
heq a b = (hsub a b) && (hsub b a)


instance Ord (HLinRel a b) where
  (<=) = hsub
  (>=) = flip hsub 

instance Eq (HLinRel a b) where
  (==) = heq

It is useful the use the direct sum of the spaces as a monoidal product.

hpar :: HLinRel a b -> HLinRel c d -> HLinRel (Either a c) (Either b d)
hpar (HLinRel mab v) (HLinRel mcd v') = HLinRel (fromBlocks [ [mab, 0], [0 , mcd]]) (vjoin [v, v']) where

hleft :: forall a b. (BEnum a, BEnum b) => HLinRel a (Either a b)
hleft = HLinRel ( i ||| (- i) ||| (konst 0 (ca,cb))) (konst 0 ca) where 
    ca = card @a
    cb = card @b  
    i = ident ca

hright :: forall a b. (BEnum a, BEnum b) => HLinRel b (Either a b)
hright = HLinRel ( i ||| (konst 0 (cb,ca)) ||| (- i) ) (konst 0 cb) where 
    ca = card @a
    cb = card @b  
    i = ident cb

htrans :: HLinRel a (Either b c) -> HLinRel (Either a b) c 
htrans (HLinRel m v) = HLinRel m v

hswap :: forall a b. (BEnum a, BEnum b) => HLinRel (Either a b) (Either b a)
hswap = HLinRel (fromBlocks [[ia ,0,0 ,-ia], [0, ib,-ib,0]]) (konst 0 (ca + cb)) where 
        ca = card @a
        cb = card @b  
        ia = ident ca
        ib = ident cb


hsum :: forall a. BEnum a => HLinRel (Either a a) a
hsum = HLinRel ( i ||| i ||| - i ) (konst 0 ca)  where 
        ca = card @a 
        i= ident ca

hdup :: forall a. BEnum a => HLinRel a (Either a a)
hdup = HLinRel (fromBlocks [[i, -i,0 ], [i, 0, -i]]) (konst 0 (ca + ca))  where 
        ca = card @a 
        i= ident ca

hdump :: HLinRel a Void
hdump = HLinRel 0 0

hlabsorb ::forall a. BEnum a => HLinRel (Either Void a) a
hlabsorb = HLinRel m v where (HLinRel m v) = hid @a 

A side note: Void causes some consternation. Void is the type with no elements and is the index type of a 0 dimensional space. It is the unit object of the monoidal product. Unfortunately by an accident of the standard Haskell definitions, actual Void is not a BEnum. So, I did a disgusting hack. Let us not discuss it more.

Circuits

Baez and Fong have an interesting paper where they describe building circuits using a categorical graphical calculus. We have the pieces to go about something similar. What we have here is a precise way in which circuit diagrams can be though of as string diagrams in a monoidal category of linear relations.

An idealized wire has two quantities associated with it, the current flowing through it and the voltage it is at.

-- a 2d space at every wire or current and voltage.
data IV = I | V deriving (Show, Enum, Bounded, Eq, Ord)

When we connect wires, the currents must be conserved and the voltages must be equal. hid and hcompose from above still achieve that. Composing two independent circuits in parallel is achieve by hpar.

Independent resistors in parallel.

We will want some basic tinker toys to work with.

A resistor in series has the same current at both ends and a voltage drop proportional to the current

resistor :: Double -> HLinRel IV IV
resistor r = HLinRel ( (2><4)  [ 1,0,-1,   0,
                                 r, 1, 0, -1]) (konst 0 2) 

Composing two resistors in parallel adds the resistance. (resistor r1) <<< (resistor r2) == resistor (r1 + r2))

A bridging resistor allows current to flow between the two branches

bridge :: Double -> HLinRel (Either IV IV) (Either IV IV)
bridge r = HLinRel (  (4><8) [ 1,0, 1,  0, -1, 0, -1,  0, -- current conservation
                               0, 1, 0, 0, 0, -1 , 0,  0, --voltage maintained left
                               0, 0, 0, 1, 0,  0,  0, -1, -- voltage maintained right
                               r, 1, 0,-1, -r,  0,  0, 0  ]) (konst 0 4)  
A bridging resistor

Composing two bridge circuits is putting the bridge resistors in parallel. The conductance G=\frac{1}{R} of resistors in parallel adds. hcompose (bridge r1) (bridge r2) == bridge 1 / (1/r1 + 1/r2).

parallel resistors compose

An open circuit allows no current to flow and ends a wire. open ~ resistor infinity

open :: HLinRel IV Void
open = HLinRel (fromList [[1,0]]) (konst 0 1)

At branching points, the voltage is maintained, but the current splits.

cmerge :: HLinRel (Either IV IV) IV
cmerge = HLinRel (fromList [[1, 0, 1, 0, -1, 0],
                           [0,1,0,0,0 ,-1  ],
                           [0,0,0,1, 0, -1]])  (konst 0 3)

This cmerge combinator could also be built using a short == bridge 0 , composing a branch with open, and then absorbing the Void away.

We can bend wires up or down by using a composition of cmerge and open.

cap :: HLinRel  (Either IV IV) Void
cap  = hcompose open cmerge

cup :: HLinRel Void (Either IV IV)
cup = hconverse cap

ground :: HLinRel IV Void
ground = HLinRel ( (1><2) [ 0 , 1 ]) (vzero 1) 

Voltage and current sources enforce current and voltage to be certain values

vsource :: Double -> HLinRel IV IV
vsource v = HLinRel ( (2><4) [ 1,0,-1,   0,
                               0, 1, 0, -1]) (fromList [0,v])  

isource :: Double -> HLinRel IV IV
isource i = HLinRel (fromList [ [1,0, -1,   0], -- current conservation
                                [1, 0, 0,  0]]) (fromList [0,i])  

Measurements of circuits proceed by probes.

type VProbe = ()
vprobe :: HLinRel IV VProbe
vprobe = HLinRel ( (2><3)  [1,0,0,
                            0,1,-1]) (konst 0 2)  

Inductors and capacitors could be included easily, but would require the entries of the HMatrix values to be polynomials in the frequency \omega, which it does not support (but it could!). We'll leave those off for another day.

We actually can determine that the rules suggested above are being followed by computation.

r20 :: HLinRel IV IV
r20 = resistor 20

main :: IO ()
main = do
            print (r20 == (hid <<< r20))
            print (r20 == r20 <<< hid)
            print (r20 == (hmeet r20 r20))
            print $ resistor 50 == r20 <<< (resistor 30)
            print $ (bridge 10) <<< (bridge 10) == (bridge 5)
            print $ v2h (h2v r20) == r20
            print $ r20 <= htop
            print $ hconverse (hconverse r20) == r20
            print $ (open <<< r20) == open

Bits and Bobbles

  • Homogenous systems are usually a bit more elegant to deal with, although a bit more unfamiliar and abstract.
  • Could make a pandas like interface for linear relations that uses numpy/scipy.sparse for the computation. All the swapping and associating is kind of fun to design, not so much to use. Labelled n-way relations are nice for users.
  • Implicit/Lazy evaluation. We should let the good solvers do the work when possible. We implemented our operations eagerly. We don't have to. By allowing hidden variables inside our relations, we can avoid the expensive linear operations until it is useful to actually compute on them.
  • Relational division = quotient spaces?
  • DSL. One of the beauties of the pointfree/categorical approach is that you avoid the need for binding forms. This makes for a very easily manipulated DSL. The transformations feel like those of ordinary algebra and you don't have to worry about the subtleties of index renaming or substitution under binders.
  • Sparse is probably really good. We have lots of identity matrices and simple rearrangements. It is very wasteful to use dense operations on these.
  • Schur complement https://en.wikipedia.org/wiki/Schur_complement are the name in the game for projecting out pieces of linear problems. We have some overlap.
  • Linear relations -> Polyhedral relations -> Convex Relations. Linear is super computable, polyhedral can blow up. Rearrange a DSL to abuse Linear programming as much as possible for queries.
  • Network circuits. There is an interesting subclass of circuits that is designed to be pretty composable.

https://en.wikipedia.org/wiki/Two-port_network Two port networks are a very useful subclass of electrical circuits. They model transmission lines fairly well, and easily composable for filter construction.

It is standard to describe these networks by giving a linear function between two variables and the other two variables. Depending on your choice of which variables depend on which, these are called the z-parameters, y-parameters, h-parameters, scattering parameters, abcd parameters. There are tables of formula for converting from one form to the others. The different parameters hold different use cases for composition and combining in parallel or series. From the perspective of linear relations this all seems rather silly. The necessity for so many descriptions and the confusing relationship between them comes from the unnecessary and overly rigid requirement of have a linear function-like relationship rather than just a general relation, which depending of the circuit may not even be available (there are degenerate configurations where two of the variables do not imply the values of the other two). A function relationship is always a lie (although a sometimes useful one), as there is always back-reaction of new connections.

-- voltage divider
divider :: Double -> Double -> HLinRel (Either IV IV) (Either IV IV)
divider r1 r2 = hcompose (bridge r2) (hpar (resistor r1) hid) 

The relation model also makes clearer how to build lumped models out of continuous ones. https://en.wikipedia.org/wiki/Lumped-element_model

https://en.wikipedia.org/wiki/Transmission_line https://en.wikipedia.org/wiki/Distributed-element_model

    null
  • Because the type indices have no connection to the actual data types (they are phantom) it is a wise idea to use smart constructors that check that the sizes of the matrices makes sense.

-- smart constructors
hLinRel :: forall a b. (BEnum a, BEnum b) => Matrix Double -> Vector Double -> Maybe (HLinRel a b) 
hLinRel m v | cols m == (ca + cb) &&  (size v == rows m)  = Just (HLinRel m v)
            |  otherwise = Nothing  where 
                 ca = card @a
                 cb = card @b  
  • Nonlinear circuits. Grobner Bases and polynomial relations?
  • Quadratic optimization under linear constraints. Can't get it to come out right yet. Clutch for Kalman filters. Nice for many formulations like least power, least action, minimum energy principles. Edit: I did more in this direction here http://www.philipzucker.com/categorical-lqr-control-with-linear-relations/
  • Quadratic Operators -> Convex operators. See last chapter of Rockafellar.
  • Duality of controllers and filters. It is well known (I think) that for ever controller algorithm there is a filter algorithm that is basically the same thing.
    • LQR - Kalman
    • Viterbi filter - Value function table
    • particle filter - Monte Carlo control
    • Extended Kalman - iLQR-ish? Use local approximation of dynamics
    • unscented kalman - ?

Refs

Failing to Bound Kissing Numbers

https://en.wikipedia.org/wiki/Kissing_number

Cody brought up the other day the kissing number problem.Kissing numbers are the number of equal sized spheres you can pack around another one in d dimensions. It’s fairly self evident that the number is 2 for 1-d and 6 for 2d but 3d isn’t so obvious and in fact puzzled great mathematicians for a while. He was musing that it was interesting that he kissing numbers for some dimensions are not currently known, despite the fact that the first order theory of the real numbers is decidable https://en.wikipedia.org/wiki/Decidability_of_first-order_theories_of_the_real_numbers

I suggested on knee jerk that Sum of Squares might be useful here. I see inequalities and polynomials and then it is the only game in town that I know anything about.

Apparently that knee jerk was not completely wrong

https://arxiv.org/pdf/math/0608426.pdf

Somehow SOS/SDP was used for bounds here. I had an impulse that the problem feels SOS-y but I do not understand their derivation.

One way the problem can be formulated is by finding or proving there is no solution to the following set of equations constraining the centers x_i of the spheres. Set the central sphere at (0,0,0,…) . Make the radii 1. Then\forall i. |x_i|^2 = 2^2 and \forall i j.  |x_i - x_j|^2 \ge 2^2

I tried a couple different things and have basically failed. I hope maybe I’ll someday have a follow up post where I do better.

So I had 1 idea on how to approach this via a convex relaxation

Make a vector x = \begin{bmatrix} x_0 & y _0 & x_1 & y _1 & x_2 & y _2 & ... \end{bmatrix} Take the outer product of this vector x^T x = X Then we can write the above equations as linear equalities and inequalities on X. If we forget that we need X to be the outer product of x (the relaxation step), this becomes a semidefinite program. Fingers crossed, maybe the solution comes back as a rank 1 matrix. Other fingers crossed, maybe the solution comes back and says it’s infeasible. In either case, we have solved our original problem.

import numpy as np
import cvxpy as cvx


d = 2
n = 6
N = d * n 
x = cvx.Variable((N+1,N+1), symmetric=True)
c = []
c += [x >> 0]
c += [x[0,0] == 1]
# x^2 + y^2 + z^2 + ... == 2^2 constraint
x1 = x[1:,1:]
for i in range(n):
    q = 0
    for j in range(d):
        q += x1[d*i + j, d*i + j]
    c += [q == 4] #[ x1[2*i + 1, 2*i + 1] + x[2*i + 2, 2*i + 2] == 4]

#c += [x1[0,0] == 2, x1[1,1] >= 0]
#c += [x1[2,2] >= 2]

# (x - x) + (y - y) >= 4
for i in range(n):    
    for k in range(i):
        q = 0
        for j in range(d):
            q += x1[d*i + j, d*i + j] + x1[d*k + j, d*k + j] - 2 * x1[d*i + j, d*k + j] # xk ^ 2 - 2 * xk * xi 
        c += [q >= 4]
print(c)
obj = cvx.Maximize(cvx.trace(np.random.rand(N+1,N+1) @ x ))
prob = cvx.Problem(obj, c)
print(prob.solve(verbose=True))
u, s, vh = np.linalg.svd(x.value)
print(s)

import matplotlib.pyplot as plt
xy = vh[0,1:].reshape(-1,2)
print(xy)
plt.scatter(xy[0], xy[1] )
plt.show()

Didn’t work though. Sigh. It’s conceivable we might do better if we start packing higher powers into x?

Ok Round 2. Let’s just ask z3 and see what it does. I’d trust z3 with my baby’s soft spot.

It solves for 5 and below. Z3 grinds to a halt on N=6 and above. It ran for days doin nothing on my desktop.

from z3 import *
import numpy as np

d = 2 # dimensions
n = 6 # number oif spheres

x = np.array([ [ Real("x_%d_%d" % (i,j))     for j in range(d) ] for i in range(n)])
print(x)

c = []
ds = np.sum(x**2, axis=1)
c += [ d2 == 4 for d2 in ds] # centers at distance 2 from origin


ds = np.sum( (x.reshape((-1,1,d)) - x.reshape((1,-1,d)))**2, axis = 2)

c += [ ds[i,j]  >= 4  for i in range(n) for j in range(i)] # spheres greater than dist 2 apart
c += [x[0,0] == 2]
print(c)
print(solve(c))

Ok. A different tact. Try to use a positivstellensatz proof. If you have a bunch of polynomial inequalities and equalities if you sum polynomial multiples of these constraints, with the inequalities having sum of square multiples, in such a way to = -1, it shows that there is no real solution to them. We have the distance from origin as equality constraint and distance from each other as an inequality constraint. I intuitively think of the positivstellensatz as deriving an impossibility from false assumptions. You can’t add a bunch of 0 and positive numbers are get a negative number, hence there is no real solution.

I have a small set of helper functions for combining sympy and cvxpy for sum of squares optimization. I keep it here along with some other cute little constructs https://github.com/philzook58/cvxpy-helpers

import cvxpy as cvx
from sympy import *
import random
'''
The idea is to use raw cvxpy and sympy as much as possible for maximum flexibility.

Construct a sum of squares polynomial using sospoly. This returns a variable dictionary mapping sympy variables to cvxpy variables.
You are free to the do polynomial operations (differentiation, integration, algerba) in pure sympy
When you want to express an equality constraint, use poly_eq(), which takes the vardict and returns a list of cvxpy constraints.
Once the problem is solved, use poly_value to get back the solution polynomials.

That some polynomial is sum of squares can be expressed as the equality with a fresh polynomial that is explicility sum of sqaures.

With the approach, we get the full unbridled power of sympy (including grobner bases!)

I prefer manually controlling the vardict to having it auto controlled by a class, just as a I prefer manually controlling my constraint sets
Classes suck.
'''


def cvxify(expr, cvxdict): # replaces sympy variables with cvx variables in sympy expr
     return lambdify(tuple(cvxdict.keys()), expr)(*cvxdict.values()) 

def sospoly(terms, name=None):
    ''' returns sum of squares polynomial using terms, and vardict mapping to cvxpy variables '''
    if name == None:
        name = str(random.getrandbits(32))
    N = len(terms)
    xn = Matrix(terms)
    Q = MatrixSymbol(name, N,N)
    p = (xn.T * Matrix(Q) * xn)[0]
    Qcvx = cvx.Variable((N,N), PSD=True)
    vardict = {Q : Qcvx} 
    return p, vardict



def polyvar(terms,name=None):
    ''' builds sumpy expression and vardict for an unknown linear combination of the terms '''
    if name == None:
        name = str(random.getrandbits(32))
    N = len(terms)
    xn = Matrix(terms)
    Q = MatrixSymbol(name, N, 1)
    p = (xn.T * Matrix(Q))[0]
    Qcvx = cvx.Variable((N,1))
    vardict = {Q : Qcvx} 
    return p, vardict

def scalarvar(name=None):
    return polyvar([1], name)

def worker(x ):
    (expr,vardict) = x
    return cvxify(expr, vardict) == 0
def poly_eq(p1, p2 , vardict):
    ''' returns a list of cvxpy constraints '''
    dp = p1 - p2
    polyvars = list(dp.free_symbols - set(vardict.keys()))
    print("hey")
    p, opt = poly_from_expr(dp, gens = polyvars, domain = polys.domains.EX) #This is brutalizing me
    print(opt)
    print("buddo")
    return [ cvxify(expr, vardict) == 0 for expr in p.coeffs()]
    '''
    import multiprocessing
    import itertools
    pool = multiprocessing.Pool()

    return pool.imap_unordered(worker, zip(p.coeffs(),  itertools.repeat(vardict)))
    '''
    
def vardict_value(vardict):
    ''' evaluate numerical values of vardict '''
    return {k : v.value for (k, v) in vardict.items()}

def poly_value(p1, vardict):
    ''' evaluate polynomial expressions with vardict'''
    return cvxify(p1, vardict_value(vardict))

if __name__ == "__main__":
    x = symbols('x')
    terms = [1, x, x**2]
    #p, cdict = polyvar(terms)
    p, cdict = sospoly(terms)
    c = poly_eq(p, (1 + x)**2 , cdict)
    print(c)
    prob = cvx.Problem(cvx.Minimize(1), c)
    prob.solve()

    print(factor(poly_value(p, cdict)))

    # global poly minimization
    vdict = {}
    t, d = polyvar([1], name='t')
    vdict.update(d)

    p, d = sospoly([1,x,x**2], name='p')
    vdict.update(d)
    constraints = poly_eq(7 + x**2 - t, p, vdict)
    obj = cvx.Maximize( cvxify(t,vdict) )
    prob = cvx.Problem(obj, constraints)
    prob.solve()
    print(poly_value(t,vdict))

and here is the attempted positivstellensatz.

import sos
import cvxpy as cvx
from sympy import *
import numpy as np

d = 2
N = 7

# a grid of a vector field. indices = (xposition, yposition, vector component)
'''xs = [ [symbols("x_%d_%d" % (i,j)) for j in range(d)] for i in range(N) ]
gens = [x for l in xs for x in l ]
xs = np.array([[poly(x,gens=gens, domain=polys.domains.EX) for x in l] for l in xs])
'''
xs = np.array([ [symbols("x_%d_%d" % (i,j)) for j in range(d)] for i in range(N) ])

c1 = np.sum( xs * xs, axis=1) - 1
c2 = np.sum((xs.reshape(-1,1,d) - xs.reshape(1,-1,d))**2 , axis=2) - 1

print(c1)
print(c2)
terms0 = [1]
terms1 = terms0 + list(xs.flatten())
terms2 = [ terms1[i]*terms1[j] for j in range(N+1) for i in range(j+1)]
#print(terms1)
#print(terms2)
vdict = {}
psatz = 0
for c in c1:
    lam, d = sos.polyvar(terms2)
    vdict.update(d)
    psatz += lam*c
for i in range(N):
    for j in range(i):
        c = c2[i,j]
        lam, d = sos.sospoly(terms2)
        vdict.update(d)
        psatz += lam*c
#print(type(psatz))
print("build constraints")
constraints = sos.poly_eq(psatz, -1, vdict)
#print("Constraints: ", len(constraints))
obj = cvx.Minimize(1) #sum([cvx.sum(v) for v in vdict.values()]))
print("build prob")
prob = cvx.Problem(obj, constraints)
print("solve")
prob.solve(verbose=True, solver= cvx.SCS)

It worked in 1-d, but did not work in 2d. At order 3 polynomials N=7, I maxed out my ram.

I also tried doing it in Julia, since sympy was killing me. Julia already has a SOS package

using JuMP
using SumOfSquares
using DynamicPolynomials
using SCS
N = 10
d = 2
@polyvar x[1:N,1:d]
X = monomials(reshape(x,d*N), 0:2)
X1 = monomials(reshape(x,d*N), 0:4)

model = SOSModel(with_optimizer(SCS.Optimizer))

acc = nothing
for t in sum(x .* x, dims=2)
    #print(t)
    p = @variable(model, [1:1], Poly(X1))
    #print(p)
    if acc != nothing
        acc += p * (t - 1)
    else
        acc = p * (t - 1)
    end
end

for i in range(1,stop=N)
    for j in range(1,stop=i-1)
        d  = x[i,:] - x[j,:]
        p = @variable(model, [1:1], SOSPoly(X))
        acc += p * (sum(d .* d) - 1)
    end
end

#print(acc)
print(typeof(acc))
@constraint(model, acc[1] == -1 )
optimize!(model)

It was faster to encode, but it’s using the same solver (SCS), so basically the same thing.

I should probably be reducing the system with respect to equality constraints since they’re already in a Groebner basis. I know that can be really important for reducing the size of your problem

I dunno.

Blah blah blah blah A bunch of unedited trash

https://github.com/peterwittek/ncpol2sdpa Peter Wittek has probably died in an avalanche? That is very sad.

These notes

https://web.stanford.edu/class/ee364b/lectures/sos_slides.pdf

Positivstullensatz.

kissing number

Review of sum of squares

minimimum sample as LP. ridiculous problem
min t
st. f(x_i) – t >= 0

dual -> one dual variable per sample point
The only dual that will be non zero is that actually selecting the minimum.

Hm. Yeah, that’s a decent analogy.

How does the dual even have a chance of knowing about poly airhtmetic?
It must be during the SOS conversion prcoess. In building the SOS constraints,
we build a finite, limittted version of polynomial multiplication
x as a matrix. x is a shift matrix.
In prpducing the characterstic polynomial, x is a shift matrix, with the last line using the polynomial
known to be zero to
eigenvectors of this matrix are zeros of the poly.

SOS does not really on polynomials persay. It relies on closure of the suqaring operaiton

maybe set one sphere just at x=0 y = 2. That breaks some symmettry

set next sphere in plane something. random plane through origin?

order y components – breaks some of permutation symmettry.

no, why not order in a random direction. That seems better for symmettry breaking

Learn Coq in Y

Edit: It’s up! https://learnxinyminutes.com/docs/coq/

I’ve been preparing a Learn X in Y tutorial for Coq. https://learnxinyminutes.com/

I’ve been telling people this and been surprised by how few people have heard of the site. It’s super quick intros to syntax and weirdness for a bunch of languages with inline code tutorials.
I think that for me, a short description of that mundane syntactic and programming constructs of coq is helpful.
Some guidance of the standard library, what is available by default. And dealing with Notation scopes, which is a pretty weird feature that most languages don’t have.
The manual actually has all this now. It’s really good. Like check this section out https://coq.inria.fr/refman/language/coq-library.html . But the manual is an intimidating documents. It starts with a BNF description of syntax and things like that. The really useful pedagogical stuff is scattered throughout it.

Anyway here is my draft (also here https://github.com/philzook58/learnxinyminutes-docs/blob/master/coq.html.markdown where the syntax highlighting isn’t so janked up). Suggestions welcome. Or if this gets accepted, you can just make pull requests


---
language: Coq
filename: learncoq.v
contributors:
    - ["Philip Zucker", "http://www.philipzucker.com/"]
---

The Coq system is a proof assistant. It is designed to build and verify mathematical proofs. The Coq system contains the functional programming language Gallina and is capable of proving properties about programs written in this language.

Coq is a dependently typed language. This means that the types of the language may depend on the values of variables. In this respect, it is similar to other related languages such as Agda, Idris, F*, Lean, and others. Via the Curry-Howard correspondence, programs, properties and proofs are formalized in the same language.

Coq is developed in OCaml and shares some syntactic and conceptual similiarity with it. Coq is a language containing many fascinating but difficult topics. This tutorial will focus on the programming aspects of Coq, rather than the proving. It may be helpful, but not necessary to learn some OCaml first, especially if you are unfamiliar with functional programming. This tutorial is based upon its OCaml equivalent

The standard usage model of Coq is to write it with interactive tool assistance, which operates like a high powered REPL. Two common such editors are the CoqIDE and Proof General Emacs mode.

Inside Proof General `Ctrl+C ` will evaluate up to your cursor.


```coq
(*** Comments ***)

(* Comments are enclosed in (* and *). It's fine to nest comments. *)

(* There are no single-line comments. *)

(*** Variables and functions ***)

(* The Coq proof assistant can be controlled and queried by a command language called 
   the vernacular. Vernacular keywords are capitalized and the commands end with a period.
   Variable and function declarations are formed with the Definition vernacular. *)

Definition x := 10.

(* Coq can sometimes infer the types of arguments, but it is common practice to annotate
   with types. *)

Definition inc_nat (x : nat) : nat := x + 1.

(* There exists a large number of vernacular commands for querying information. 
   These can be very useful. *)

Compute (1 + 1). (* 2 : nat *) (* Compute a result. *)

Check tt. (* tt : unit *) (* Check the type of an expressions *)

About plus. (* Prints information about an object *)

(* Print information including the definition *)
Print true. (* Inductive bool : Set := true : Bool | false : Bool *)  

Search nat. (* Returns a large list of nat related values *)
Search "_ + _". (* You can also search on patterns *)
Search (?a -> ?a -> bool). (* Patterns can have named parameters  *)
Search (?a * ?a).

(* Locate tells you where notation is coming from. Very helpful when you encounter
   new notation. *)
Locate "+". 

(* Calling a function with insufficient number of arguments
   does not cause an error, it produces a new function. *)
Definition make_inc x y := x + y. (* make_inc is int -> int -> int *)
Definition inc_2 := make_inc 2.   (* inc_2 is int -> int *)
Compute inc_2 3. (* Evaluates to 5 *)

(* Definitions can be chained with "let ... in" construct.
   This is roughly the same to assigning values to multiple
   variables before using them in expressions in imperative
   languages. *)
Definition add_xy : nat := let x := 10 in
                         let y := 20 in
                         x + y.


(* Pattern matching is somewhat similar to switch statement in imperative
   languages, but offers a lot more expressive power. *)
Definition is_zero (x : nat) :=
    match x with
    | 0 => true
    | _ => false  (* The "_" pattern means "anything else". *)
    end.


(* You can define recursive function definition using the Fixpoint vernacular.*)
Fixpoint factorial n := match n with
                        | 0 => 1
                        | (S n') => n * factorial n'
                        end.


(* Function application usually doesn't need parentheses around arguments *)
Compute factorial 5. (* 120 : nat *)

(* ...unless the argument is an expression. *)
Compute factorial (5-1). (* 24 : nat *)

(* You can define mutually recursive functions using "with" *)
Fixpoint is_even (n : nat) : bool := match n with
  | 0 => true
  | (S n) => is_odd n
end with
  is_odd n := match n with
  | 0 => false
  | (S n) => is_even n
              end.

(* As Coq is a total programming language, it will only accept programs when it can
   understand they terminate. It can be most easily seen when the recursive call is
   on a pattern matched out subpiece of the input, as then the input is always decreasing
   in size. Getting Coq to understand that functions terminate is not always easy. See the
   references at the end of the artice for more on this topic. *)

(* Anonymous functions use the following syntax: *)

Definition my_square : nat -> nat := fun x => x * x.

Definition my_id (A : Type) (x : A) : A := x.
Definition my_id2 : forall A : Type, A -> A := fun A x => x.
Compute my_id nat 3. (* 3 : nat *)

(* You can ask Coq to infer terms with an underscore *)
Compute my_id _ 3. 

(* An implicit argument of a function is an argument which can be inferred from contextual
   knowledge. Parameters enclosed in {} are implicit by default *)

Definition my_id3 {A : Type} (x : A) : A := x.
Compute my_id3 3. (* 3 : nat *)

(* Sometimes it may be necessary to turn this off. You can make all arguments explicit
   again with @ *)
Compute @my_id3 nat 3.

(* Or give arguments by name *)
Compute my_id3 (A:=nat) 3.

(*** Notation ***)

(* Coq has a very powerful Notation system that can be used to write expressions in more
   natural forms. *)
Compute Nat.add 3 4. (* 7 : nat *)
Compute 3 + 4. (* 7 : nat *)

(* Notation is a syntactic transformation applied to the text of the program before being
   evaluated. Notation is organized into notation scopes. Using different notation scopes allows for a weak notion of overloading. *)

(* Imports the Zarith module containing definitions related to the integers Z *)
Require Import ZArith. 

(* Notation scopes can be opened *)
Open Scope Z_scope.

(* Now numerals and addition are defined on the integers. *)
Compute 1 + 7. (* 8 : Z *)

(* Integer equality checking *)
Compute 1 =? 2. (* false : bool *) 

(* Locate is useful for finding the origin and definition of notations *)
Locate "_ =? _". (* Z.eqb x y : Z_scope *) 
Close Scope Z_scope.

(* We're back to nat being the default interpetation of "+" *)
Compute 1 + 7. (* 8 : nat *)

(* Scopes can also be opened inline with the shorthand % *)
Compute (3 * -7)%Z. (* -21%Z : Z *)

(* Coq declares by default the following interpretation scopes: core_scope, type_scope, 
   function_scope, nat_scope, bool_scope, list_scope, int_scope, uint_scope. You may also
   want the numerical scopes Z_scope (integers) and Q_scope (fractions) held in the ZArith
   and QArith module respectively. *)

(* You can print the contents of scopes *)
Print Scope nat_scope.
(*
Scope nat_scope
Delimiting key is nat
Bound to classes nat Nat.t
"x 'mod' y" := Nat.modulo x y
"x ^ y" := Nat.pow x y
"x ?= y" := Nat.compare x y
"x >= y" := ge x y
"x > y" := gt x y
"x =? y" := Nat.eqb x y
"x  a
                                                  end.

(* A destructuring let is available if a pattern match is irrefutable *)
Definition my_fst2 {A B : Type} (x : A * B) : A := let (a,b) := x in
                                                   a.

(*** Lists ***)

(* Lists are built by using cons and nil or by using notation available in list_scope. *)
Compute cons 1 (cons 2 (cons 3 nil)). (*  (1 :: 2 :: 3 :: nil)%list : list nat *)
Compute (1 :: 2 :: 3 :: nil)%list. 

(* There is also list notation available in the ListNotations modules *)
Require Import List.
Import ListNotations. 
Compute [1 ; 2 ; 3]. (* [1; 2; 3] : list nat *)


(* 
There are a large number of list manipulation functions available, lncluding:

• length
• head : first element (with default) 
• tail : all but first element
• app : appending
• rev : reverse
• nth : accessing n-th element (with default)
• map : applying a function
• flat_map : applying a function returning lists 
• fold_left : iterator (from head to tail)
• fold_right : iterator (from tail to head) 

 *)

Definition my_list : list nat := [47; 18; 34].

Compute List.length my_list. (* 3 : nat *)
(* All functions in coq must be total, so indexing requires a default value *)
Compute List.nth 1 my_list 0. (* 18 : nat *) 
Compute List.map (fun x => x * 2) my_list. (* [94; 36; 68] : list nat *)
Compute List.filter (fun x => Nat.eqb (Nat.modulo x 2) 0) my_list. (*  [18; 34] : list nat *)
Compute (my_list ++ my_list)%list. (*  [47; 18; 34; 47; 18; 34] : list nat *)

(*** Strings ***)

Require Import Strings.String.

Open Scope string_scope.

(* Use double quotes for string literals. *)
Compute "hi"%string.

(* Strings can be concatenated with the "++" operator. *)
Compute String.append "Hello " "World". (* "Hello World" : string *)
Compute "Hello " ++ "World". (* "Hello World" : string *)

(* Strings can be compared for equality *)
Compute String.eqb "Coq is fun!"%string "Coq is fun!"%string. (* true : bool *)
Compute ("no" =? "way")%string. (* false : bool *)

Close Scope string_scope.

(*** Other Modules ***)

(* Other Modules in the standard library that may be of interest:

• Logic : Classical logic and dependent equality
• Arith : Basic Peano arithmetic
• PArith : Basic positive integer arithmetic
• NArith : Basic binary natural number arithmetic
• ZArith : Basic relative integer arithmetic
• Numbers : Various approaches to natural, integer and cyclic numbers (currently axiomatically and on top of 2^31 binary words)
• Bool : Booleans (basic functions and results)
• Lists : Monomorphic and polymorphic lists (basic functions and results), Streams (infinite sequences
defined with co-inductive types)
• Sets : Sets (classical, constructive, finite, infinite, power set, etc.)
• FSets : Specification and implementations of finite sets and finite maps (by lists and by AVL trees)
• Reals : Axiomatization of real numbers (classical, basic functions, integer part, fractional part, limit, derivative, Cauchy series, power series and results,...)
• Relations : Relations (definitions and basic results)
• Sorting : Sorted list (basic definitions and heapsort correctness)
• Strings : 8-bits characters and strings
• Wellfounded : Well-founded relations (basic results)
 *)

(*** User-defined data types ***)

(* Because Coq is dependently typed, defining type aliases is no different than defining
   an alias for a value. *)

Definition my_three : nat := 3.
Definition my_nat : Type := nat.

(* More interesting types can be defined using the Inductive vernacular. Simple enumeration
   can be defined like so *)
Inductive ml := OCaml | StandardML | Coq.
Definition lang := Coq.  (* Has type "ml". *)

(* For more complicated types, you will need to specify the types of the constructors. *)

(* Type constructors don't need to be empty. *)
Inductive my_number := plus_infinity
                     | nat_value : nat -> my_number.
Compute nat_value 3. (* nat_value 3 : my_number *)


(* Record syntax is sugar for tuple-like types. It defines named accessor functions for
   the components *)
Record Point2d (A : Set) := mkPoint2d { x2 : A ; y2 : A }. 
Definition mypoint : Point2d nat :=  {| x2 := 2 ; y2 := 3 |}.
Compute x2 nat mypoint. (* 2 : nat *)
Compute mypoint.(x2 nat). (* 2 : nat *) 

(* Types can be parameterized, like in this type for "list of lists
   of anything". 'a can be substituted with any type. *)
Definition list_of_lists a := list (list a).
Definition list_list_nat := list_of_lists nat.

(* Types can also be recursive. Like in this type analogous to
   built-in list of naturals. *)

Inductive my_nat_list := EmptyList | NatList : nat -> my_nat_list -> my_nat_list.
Compute NatList 1 EmptyList. (*  NatList 1 EmptyList : my_nat_list *)

(** Matching type constructors **)

Inductive animal := Dog : string -> animal | Cat : string -> animal.

Definition say x :=
    match x with
    | Dog x => (x ++ " says woof")%string
    | Cat x => (x ++ " says meow")%string
    end.

Compute say (Cat "Fluffy"). (* "Fluffy says meow". *)

(** Traversing data structures with pattern matching **)

(* Recursive types can be traversed with pattern matching easily.
   Let's see how we can traverse a data structure of the built-in list type.
   Even though the built-in cons ("::") looks like an infix operator,
   it's actually a type constructor and can be matched like any other. *)
Fixpoint sum_list l :=
    match l with
    | [] => 0
    | head :: tail => head + (sum_list tail)
    end.

Compute sum_list [1; 2; 3]. (* Evaluates to 6 *)


(*** A Taste of Proving ***)
(* Explaining the proof language is out of scope for this tutorial, but here is a taste to
   whet your appetite. Check the resources below for more. *)

(* A fascinating feature of dependently type based theorem provers is that the same
  primitive constructs underly the proof language as the programming features.
  For example, we can write and prove the proposition A and B implies A in raw Gallina *)

Definition my_theorem : forall A B, A /\ B -> A := fun A B ab => match ab with
                                            | (conj a b) => a
                                                       end.

(* Or we can prove it using tactics. Tactics are a macro language to help build proof terms
   in a more natural style and automate away some drudgery. *)
Theorem my_theorem2 : forall A B, A /\ B -> A.
Proof.
  intros A B ab.  destruct ab as [ a b ]. apply a.
Qed.

(* We can prove easily prove simple polynomial equalities using the automated tactic ring. *)
Require Import Ring.
Require Import Arith.
Theorem simple_poly : forall (x : nat), (x + 1) * (x + 2) = x * x + 3 * x + 2.
  Proof. intros. ring. Qed.

(* Here we prove the closed form for the sum of all numbers 1 to n using induction *)

Fixpoint sumn (n : nat) : nat :=
  match n with
  | 0 => 0
  | (S n') => n + (sumn n')
  end.

Theorem sum_formula : forall n, 2 * (sumn n) = (n + 1) * n.
Proof. intros. induction n.
       - reflexivity. (* 0 = 0 base case *)
       - simpl. ring [IHn]. (* induction step *)
Qed.
```
                                
With this we have only scratched the surface of Coq. It is a massive ecosystem with many interesting and peculiar topics leading all the way up to modern research.

## Further reading

* [The Coq reference manual](https://coq.inria.fr/refman/)
* [Software Foundations](https://softwarefoundations.cis.upenn.edu/)
* [Certfied Programming with Dependent Types](http://adam.chlipala.net/cpdt/)
* [Mathematical Components](https://math-comp.github.io/mcb/)
* [Coq'Art: The Calculus of Inductive Constructions](http://www.cse.chalmers.se/research/group/logic/TypesSS05/resources/coq/CoqArt/)
* [FRAP](http://adam.chlipala.net/frap/)


Bonus. An uneditted list of tactics. You’d probably prefer https://pjreddie.com/coq-tactics/



(*** Tactics ***)
 (* Although we won't explain their use in detail, here is a list of common tactics. *)

(* 

    * exact
    * simpl
    * intros
    * apply
    * assert
    * destruct
    * induction
    * reflexivity
    * rewrite
    * inversion
    * injection
    * discriminate
    * fold
    * unfold


        Tacticals
    * try
    * ;
        * repeat
        *
      

    
 Automatic
    * auto
    * eauto
    * tauto
    * ring
    * ring_simplify
    * psatz
    * lia
    * ria


LTac is a logic programming scripting language for tactics
     
       From Tatics chapter of LF
        intros: move hypotheses/variables from goal to context
reflexivity: finish the proof (when the goal looks like e = e)
apply: prove goal using a hypothesis, lemma, or constructor
apply... in H: apply a hypothesis, lemma, or constructor to a hypothesis in the context (forward reasoning)
apply... with...: explicitly specify values for variables that cannot be determined by pattern matching
simpl: simplify computations in the goal
simpl in H: ... or a hypothesis
rewrite: use an equality hypothesis (or lemma) to rewrite the goal
rewrite ... in H: ... or a hypothesis
symmetry: changes a goal of the form t=u into u=t
symmetry in H: changes a hypothesis of the form t=u into u=t
unfold: replace a defined constant by its right-hand side in the goal
unfold... in H: ... or a hypothesis
destruct... as...: case analysis on values of inductively defined types
destruct... eqn:...: specify the name of an equation to be added to the context, recording the result of the case analysis
induction... as...: induction on values of inductively defined types
injection: reason by injectivity on equalities between values of inductively defined types
discriminate: reason by disjointness of constructors on equalities between values of inductively defined types
assert (H: e) (or assert (e) as H): introduce a "local lemma" e and call it H
                                              generalize dependent x: move the variable x (and anything else that depends on it) from the context back to an explicit hypothesis in the goal formula
                                                                                                                                                                                        
clear H: Delete hypothesis H from the context.
subst x: For a variable x, find an assumption x = e or e = x in the context, replace x with e throughout the context and current goal, and clear the assumption.
subst: Substitute away all assumptions of the form x = e or e = x (where x is a variable).
rename... into...: Change the name of a hypothesis in the proof context. For example, if the context includes a variable named x, then rename x into y will change all occurrences of x to y.
assumption: Try to find a hypothesis H in the context that exactly matches the goal; if one is found, behave like apply H.
contradiction: Try to find a hypothesis H in the current context that is logically equivalent to False. If one is found, solve the goal.
constructor: Try to find a constructor c (from some Inductive definition in the current environment) that can be applied to solve the current goal. If one is found, behave like apply c.
                                                                                       (* Dependent types. Using dependent types for programming tasks tends to be rather unergonomic in Coq. 
We briefly mention here as an advanced topic that there exists a more sophistictaed match statement that is needed for dependently typed. See for example the "Convoy" pattern.
*)

(*** Other topics ***)

(* Dependently Typed Programming - Most of the above syntax has its equivalents in OCaml. Coq also has the capability for full dependently typed programming. There is an extended pattern matching syntax available for this purpose.
 
   Extraction - Coq can be extracted to OCaml and Haskell code for their more performant runtimes and ecosystems
   Modules / TypeClasses - Modules and Typeclasses are methods for organizing code. They allow a different form of overloading than Notation
   Setoids - 
   Termination - Gallina is a total functional programming language. It will not allow you to write functions that do not obviously terminate. For functions that do terminate but non-obviously, it requires some work to get Coq to understand this.
   Coinductive - Coinductive types such as streams are possibly infinite values that stay productive.


 *)
            

Neural Networks with Weighty Lenses (DiOptics?)

I wrote a while back how you can make a pretty nice DSL for reverse mode differentiation based on the same type as Lens. I’d heard some interesting rumblings on the internet around these ideas and so was revisiting them.

type Lens s t a b = s -> (a, b -> t)
type AD x dx y dy = x -> (y, dy -> dx)

Composition is defined identically for reverse mode just as it is for lens.

The forward computation shares info with the backwards differential propagation, which corresponds to a transposed Jacobian

After chewing on it a while, I realized this really isn’t that exotic. How it works is that you store the reverse mode computation graph, and all necessary saved data from the forward pass in the closure of the (dy -> dx). I also have a suspicion that if you defunctionalized this construction, you’d get the Wengert tape formulation of reverse mode ad.

Second, Lens is just a nice structure for bidirectional computation, with one forward pass and one backward pass which may or may not be getting/setting. There are other examples for using it like this.

It is also pretty similar to the standard “dual number” form type FAD x dx y dy = (x,dx)->(y,dy) for forward mode AD. We can bring the two closer by a CPS/Yoneda transformation and then some rearrangement.

     x -> (y, dy -> dx) 
==>  x -> (y, forall s. (dx -> s) -> (dy -> s))
==>  forall s. (x, dx -> s) -> (y, dx -> s) 

and meet it in the middle with

(x,dx) -> (y,dy)
==> forall s. (x, s -> dx) -> (y, s -> dy)

I ended the previous post somewhat unsatisfied by how ungainly writing that neural network example was, and I called for Conal Elliot’s compiling to categories plugin as a possible solution. The trouble is piping the weights all over the place. This piping is very frustrating in point-free form, especially when you know it’d be so trivial pointful. While the inputs and outputs of layers of the network compose nicely (you no longer need to know about the internal computations), the weights do not. As we get more and more layers, we get more and more weights. The weights are in some sense not as compositional as the inputs and outputs of the layers, or compose in a different way that you need to maintain access to.

I thought of a very slight conceptual twist that may help.

The idea is we keep the weights out to the side in their own little type parameter slots. Then we define composition such that it composes input/outputs while tupling the weights. Basically we throw the repetitive complexity appearing in piping the weights around into the definition of composition itself.

These operations are easily seen as 2 dimensional diagrams.

Three layers composed, exposing the weights from all layers
The 2-D arrow things can be built out of the 1-d arrows of the original basic AD lens by bending the weights up and down. Ultimately they are describing the same thing

Here’s the core reverse lens ad combinators

import Control.Arrow ((***))

type Lens'' a b = a -> (b, b -> a)

comp :: (b -> (c, (c -> b))) -> (a -> (b, (b -> a))) -> (a -> (c, (c -> a)))
comp f g x = let (b, dg) = g x in
             let (c, df) = f b in
             (c, dg . df)

id' :: Lens'' a a
id' x = (x, id) 

relu' :: (Ord a, Num a) => Lens'' a a
relu' = \x -> (frelu x, brelu x) where
        frelu x | x > 0 = x
                | otherwise = 0
        brelu x dy | x > 0 = dy
                   | otherwise = 0

add' :: Num a => Lens'' (a,a) a 
add' = \(x,y) -> (x + y, \ds -> (ds, ds))

dup' :: Num a => Lens'' a (a,a)
dup' = \x -> ((x,x), \(dx,dy) -> dx + dy)

sub' :: Num a => Lens'' (a,a) a 
sub' = \(x,y) -> (x - y, \ds -> (ds, -ds))

mul' :: Num a => Lens'' (a,a) a 
mul' = \(x,y) -> (x * y, \dz -> (dz * y, x * dz))

recip' :: Fractional a => Lens'' a a 
recip' = \x-> (recip x, \ds -> - ds / (x * x))

div' :: Fractional a => Lens'' (a,a) a 
div' = (\(x,y) -> (x / y, \d -> (d/y,-x*d/(y * y))))

sin' :: Floating a => Lens'' a a
sin' = \x -> (sin x, \dx -> dx * (cos x))

cos' :: Floating a => Lens'' a a
cos' = \x -> (cos x, \dx -> -dx * (sin x))

pow' :: Num a => Integer -> Lens'' a a
pow' n = \x -> (x ^ n, \dx -> (fromInteger n) * dx * x ^ (n-1)) 

--cmul :: Num a => a -> Lens' a a
--cmul c = lens (* c) (\x -> \dx -> c * dx)

exp' :: Floating a => Lens'' a a
exp' = \x -> let ex = exp x in
                      (ex, \dx -> dx * ex)

fst' :: Num b => Lens'' (a,b) a
fst' = (\(a,b) -> (a, \ds -> (ds, 0)))

snd' :: Num a => Lens'' (a,b) b
snd' = (\(a,b) -> (b, \ds -> (0, ds)))

-- some monoidal combinators
swap' :: Lens'' (a,b) (b,a)
swap' = (\(a,b) -> ((b,a), \(db,da) -> (da, db)))

assoc' :: Lens'' ((a,b),c) (a,(b,c))
assoc' = \((a,b),c) -> ((a,(b,c)), \(da,(db,dc)) -> ((da,db),dc))

assoc'' :: Lens'' (a,(b,c)) ((a,b),c)
assoc'' = \(a,(b,c)) -> (((a,b),c), \((da,db),dc)->  (da,(db,dc)))

par' :: Lens'' a b -> Lens'' c d -> Lens'' (a,c) (b,d)
par' l1 l2 = l3 where
    l3 (a,c) = let (b , j1) = l1 a in
               let (d, j2) = l2 c in
               ((b,d) , j1 *** j2) 
first' :: Lens'' a b -> Lens'' (a, c) (b, c)
first' l = par' l id'

second' :: Lens'' a b -> Lens'' (c, a) (c, b)
second' l = par' id' l

labsorb :: Lens'' ((),a) a
labsorb (_,a) = (a, \a' -> ((),a'))

labsorb' :: Lens'' a ((),a)
labsorb' a = (((),a), \(_,a') -> a')

rabsorb :: Lens'' (a,()) a
rabsorb = comp labsorb swap'

And here are the two dimensional combinators. I tried to write them point-free in terms of the combinators above to demonstrate that there is no monkey business going on. We

type WAD' w w' a b = Lens'' (w,a) (w',b)
type WAD'' w a b = WAD' w () a b -- terminate the weights for a closed network
{- For any monoidal category we can construct this composition? -}
-- horizontal composition
hcompose :: forall w w' w'' w''' a b c. WAD' w' w'' b c -> WAD' w w''' a b -> WAD' (w',w) (w'',w''') a c
hcompose f g = comp f' g' where 
               f' :: Lens'' ((w',r),b) ((w'',r),c)
               f' = (first' swap') `comp` assoc'' `comp` (par' id' f) `comp` assoc' `comp`  (first' swap') 
               g' :: Lens'' ((r,w),a) ((r,w'''),b)
               g' = assoc'' `comp` (par' id' g) `comp` assoc' 



rotate :: WAD' w w' a b -> WAD' a b w w'                                      
rotate f = swap' `comp` f `comp` swap'

-- vertical composition of weights
vcompose :: WAD' w'  w'' c d -> WAD' w w' a b -> WAD' w w'' (c, a) (d, b)
vcompose f g = rotate (hcompose (rotate f)  (rotate g) )                             

-- a double par.
diagpar :: forall w w' a b w'' w''' c d. WAD' w  w' a b -> WAD' w'' w''' c d 
           -> WAD' (w,w'') (w',w''') (a, c) (b, d)
diagpar f g = t' `comp` (par' f g) `comp` t where
                t :: Lens'' ((w,w''),(a,c)) ((w,a), (w'',c)) -- yikes. just rearrangements.
                t =  assoc'' `comp` (second' ((second' swap') `comp` assoc' `comp` swap')) `comp` assoc'
                t' :: Lens'' ((w',b), (w''',d)) ((w',w'''),(b,d)) -- the tranpose of t
                t' =  assoc'' `comp` (second'  ( swap'  `comp` assoc'' `comp` (second' swap')))  `comp` assoc'

id''' :: WAD' () () a a
id''' = id'






-- rotate:: WAD' w a a w
-- rotate = swap'

liftIO :: Lens'' a b -> WAD' w w a b
liftIO = second'

liftW :: Lens'' w w' -> WAD' w w' a a
liftW = first'


wassoc' = liftW assoc' 
wassoc'' = liftW assoc'' 

labsorb'' :: WAD' ((),w) w a a
labsorb'' = first' labsorb

labsorb''' :: WAD' w ((),w) a a
labsorb''' = first' labsorb'

wswap' :: WAD' (w,w') (w',w) a a
wswap' = first' swap'
-- and so on we can lift all combinators

I wonder if this is actually nice?

I asked around and it seems like this idea may be what davidad is talking about when he refers to dioptics

http://events.cs.bham.ac.uk/syco/strings3-syco5/slides/dalrymple.pdf

Perhaps this will initiate a convo.

Edit: He confirms that what I’m doing appears to be a dioptic. Also he gave a better link http://events.cs.bham.ac.uk/syco/strings3-syco5/papers/dalrymple.pdf

He is up to some interesting diagrams

https://twitter.com/davidad/status/1179760373030801408?s=20

Bits and Bobbles

  • Does this actually work or help make things any better?
  • Recurrent neural nets flip my intended role of weights and inputs.
  • Do conv-nets naturally require higher dimensional diagrams?
  • This weighty style seems like a good fit for my gauss seidel and iterative LQR solvers. A big problem I hit there was getting all the information to the outside, which is a similar issue to getting the weights around in a neural net.

Gröbner Bases and Optics

Geometrical optics is a pretty interesting topic. It really is almost pure geometry/math rather than physics.

Huygens principle says that we can compute the propagation of a wave by considering the wavelets produced by each point on the wavefront separately.

In physical optics, this corresponds to the linear superposition of the waves produced at each point by a propagator function \int dx' G(x,x'). In geometrical optics, which was Huygens original intent I think (these old school guys were VERY geometrical), this is the curious operation of taking the geometrical envelope of the little waves produced by each point.

The gist of Huygens principles. Ripped from wikipedia

https://en.wikipedia.org/wiki/Envelope_(mathematics) The envelope is an operation on a family of curves. You can approximate it by a finite difference procedure. Take two subsequent curves close together in the family, find their intersection. Keep doing that and the join up all the intersections. This is roughly the approach I took in this post http://www.philipzucker.com/elm-eikonal-sol-lewitt/

Taking the envelope of a family of lines. Ripped from wikipedia

You can describe a geometrical wavefront implicitly with an equations \phi(x,y) = 0. Maybe the wavefront is a circle, or a line, or some wacky shape.

The wavelet produced by the point x,y after a time t is described implicitly by d(\vec{x},\vec{x'})^2 - t^2 = (x-x')^2 + (y-y')^2 - t^2 = 0.

This described a family of curves, the circles produced by the different points of the original wavefront. If you take the envelope of this family you get the new wavefront at time t.

How do we do this? Grobner bases is one way if we make \phi a polynomial equation. For today’s purposes, Grobner bases are a method for solving multivariate polynomial equations. Kind of surprising that such a thing even exists. It’s actually a guaranteed terminating algorithm with horrific asymptotic complexity. Sympy has an implementation. For more on Grobner bases, the links here are useful http://www.philipzucker.com/dump-of-nonlinear-algebra-algebraic-geometry-notes-good-links-though/. Especially check out the Cox Little O’Shea books

The algorithm churns on a set of multivariate polynomials and spits out a new set that is equivalent in the sense that the new set is equal to zero if and only if the original set was. However, now (if you ask for the appropriate term ordering) the polynomials are organized in such a way that they have an increasing number of variables in them. So you solve the 1-variable equation (easy), and substitute into the 2 variable equation. Then that is a 1-variable equation, which you solve (easy) and then you substitute into the three variable equation, and so on. It’s analogous to gaussian elimination.

So check this out

from sympy import *


x1, y1, x2, y2, dx, dy = symbols('x1, y1, x2, y2, dx, dy')

def dist(x,y,d):
    return x**2 + y**2 - d**2

e1 = dist(x1,y1,2) # the original circle of radius 2
e2 = dist(x1-x2 ,y1 - y2 , 1) # the parametrized wavefront after time 1


# The two envelope conditions.
e3 = diff(e1,x1)*dx + diff(e1,y1)*1
e4 = diff(e2,x1)*dx + diff(e2,y1)*1


envelope = groebner([e1,e2,e3,e4], y1, x1, dx, dy, x2, y2, order='lex')[-1]
plot_implicit(e1, show=False)
plot_implicit(envelope, show = True)

The envelope conditions can be found by introducing two new differential variables dx, and dy. They are constrained to lie tangent to the point on the original circle by the differential equation e3, and then we require that the two subsequent members of the curve family intersect by the equation e4. Hence we get the envelope. Ask for the Grobner basis with that variable ordering gives us an implicit equations for x2, y2 with no mention of the rest if we just look at the last equation of the Grobner basis.

I set arbitrarily dy = 1 because the overall scale of them does not matter, only the direction. If you don’t do this, the final equation is scaled homogenously in dy.

This does indeed show the two new wavefronts at radius 1 and radius 3.

Original circle radius = 2

x1**2 + y1**2 – 4 = 0
Evolved circles found via grobner basis.


(x2**2 + y2**2 – 9)*(x2**2 + y2**2 – 1) = 0

Here’s a different one of a parabola using e1 =  y1 – x1 + x1**2

Original curve y1 – x1 + x1**2 = 0
After 1 time step.




16*x2**6 – 48*x2**5 + 16*x2**4*y2**2 + 32*x2**4*y2 + 4*x2**4 – 32*x2**3*y2**2 – 64*x2**3*y2 + 72*x2**3 + 32*x2**2*y2**3 + 48*x2**2*y2 – 40*x2**2 – 32*x2*y2**3 + 16*x2*y2**2 – 16*x2*y2 – 4*x2 + 16*y2**4 + 32*y2**3 – 20*y2**2 – 36*y2 – 11 = 0

The weird lumpiness here is plot_implicit’s inability to cope, not the actually curve shape Those funky prongs are from a singularity that forms as the wavefront folds over itself.

I tried using a cubic curve x**3 and the grobner basis algorithm seems to crash my computer. 🙁 Perhaps this is unsurprising. https://en.wikipedia.org/wiki/Elliptic_curve ?

I don’t know how to get the wavefront to go in only 1 direction? As is, it is propagating into the past and the future. Would this require inequalities? Sum of squares optimization perhaps?

Edit:

It’s been suggested on reddit that I’d have better luck using other packages, like Macaulay2, MAGMA, or Singular. Good point

Also it was suggested I use the Dixon resultant, for which there is an implementation in sympy, albeit hidden. Something to investigate

https://github.com/sympy/sympy/blob/master/sympy/polys/multivariate_resultants.py

https://nikoleta-v3.github.io/blog/2018/06/05/resultant-theory.html

Another interesting angle might be to try to go numerical with a homotopy continuation method with phcpy

http://homepages.math.uic.edu/~jan/phcpy_doc_html/welcome.html

https://www.semion.io/doc/solving-polynomial-systems-with-phcpy

or pybertini https://ofloveandhate-pybertini.readthedocs.io/en/feature-readthedocs_integration/intro.html

Concolic Weakest Precondition is Kind of Like a Lens

That’s a mouthful.

Lens are described as functional getters and setters. The simple lens type is

type Lens a b = a -> (b, b -> a)

. The setter is

a->b

and the getter is

a -> b -> a

This type does not constrain lenses to obey the usual laws of getters and setters. So we can use/abuse lens structures for nontrivial computations that have forward and backwards passes that share information. Jules Hedges is particular seems to be a proponent for this idea.

I’ve described before how to encode reverse mode automatic differentiation in this style. I have suspicions that you can make iterative LQR and guass-seidel iteration have this flavor too, but I’m not super sure. My attempts ended somewhat unsatisfactorily a whiles back but I think it’s not hopeless. The trouble was that you usually want the whole vector back, not just its ends.

I’ve got another example in imperative program analysis that kind of makes sense and might be useful though. Toy repo here: https://github.com/philzook58/wp-lens

In program analysis it sometimes helps to run a program both concretely and symbolically. Concolic = CONCrete / symbOLIC. Symbolic stuff can slowly find hard things and concrete execution just sprays super fast and can find the dumb things really quick.  

We can use a lens structure to organize a DSL for describing a simple imperative language

The forward pass is for the concrete execution. The backward pass is for transforming the post condition to a pre condition in a weakest precondition analysis. Weakest precondition semantics is a way of specifying what is occurring in an imperative language. It tells how each statement transforms post conditions (predicates about the state after the execution) into pre conditions (predicates about before the execution).  The concrete execution helps unroll loops and avoid branching if-then-else behavior that would make the symbolic stuff harder to process. I’ve been flipping through Djikstra’s book on this. Interesting stuff, interesting man.

I often think of a state machine as a function taking s -> s. However, this is kind of restrictive. It is possible to have heterogenous transformations s -> s’. Why not? I think I am often thinking about finite state machines, which we really don’t intend to have a changing state size. Perhaps we allocated new memory or something or brought something into or out of scope. We could model this by assuming the memory was always there, but it seems wasteful and perhaps confusing. We need to a priori know everything we will need, which seems like it might break compositionally.

We could model our language making some data type like
data Imp = Skip | Print String | Assign String Expr | Seq Imp Imp | ...
and then build an interpreter

interp :: Imp -> s -> s'

Imp.

But we can also cut out the middle man and directly define our language using combinators.

type Stmt s s' = s ->s'

To me this has some flavor of a finally tagless style.


Likewise for expressions. Expressions evaluate to something in the context of the state (they can lookup variables), so let’s just use

type Expr s a = s -> a

And, confusingly (sorry), I think it makes sense to use Lens in their original getter/setter intent for variables. So Lens structure is playing double duty.

type Var s a = Lens' s a

With that said, here we go.


type Stmt s s' = s -> s' 
type Lens' a b = a -> (b, b -> a)
set l s a = let (_, f) = l s in f a

type Expr s a = s -> a
type Var s a = Lens' s a

skip :: Stmt s s
skip = id

sequence :: Stmt s s' -> Stmt s' s'' -> Stmt s s''
sequence = flip (.)

assign :: Var s a -> Expr s a -> Stmt s s
assign v e = \s -> set v s (e s)

(===) :: Var s a -> Expr s a -> Stmt s s
v === e = assign v e

ite :: Expr s Bool -> Stmt s s' -> Stmt s s' -> Stmt s s'
ite e stmt1 stmt2 = \s -> if (e s) then stmt1 s else stmt2 s

while :: Expr s Bool -> Stmt s s -> Stmt s s
while e stmt = \s -> if (e s) then ((while e stmt) (stmt s)) else s

assert :: Expr s Bool -> Stmt s s  
assert e = \s -> if (e s) then s else undefined 

abort :: Stmt s s'  
abort = const undefined

Weakest precondition can be done similarly, instead we start from the end and work backwards

Predicates are roughly sets. A simple type for sets is

type Pred s = s -> Bool 
Now, this doesn’t have much deductive power, but I think it demonstrates the principles simply. We could replace Pred with perhaps an SMT solver expression, or some data type for predicates, for which we’ll need to implement things like substitution. Let’s not today.

A function

a -> b 
is equivalent to
forall c. (b -> c) -> (a -> c)
. This is some kind of CPS / Yoneda transformation thing. A state transformer
s -> s'
to predicate transformer
(s' -> Bool) -> (s -> Bool)
is somewhat evocative of that. I’m not being very precise here at all.

Without further ado, here’s how I think a weakest precondition looks roughly.


type Lens' a b = a -> (b, b -> a)
set l s a = let (_, f) = l s in f a

type Expr s a = s -> a
type Var s a = Lens' s a
type Pred s = s -> Bool
type Stmt s s' = Pred s' -> Pred s 

skip :: Stmt s s
skip = \post -> let pre = post in pre -- if

sequence :: Stmt s s' -> Stmt s' s'' -> Stmt s s''
sequence = (.)

assign :: Var s a -> Expr s a -> Stmt s s
assign v e = \post -> let pre s = post (set v s (e s)) in pre

(===) :: Var s a -> Expr s a -> Stmt s s
v === e = assign v e

ite :: Expr s Bool -> Stmt s s' -> Stmt s s' -> Stmt s s'
ite e stmt1 stmt2 = \post -> let pre s = if (e s) then (stmt1 post) s else (stmt2 post) s in pre

abort :: Stmt s s'  
abort = \post -> const False

assert :: Expr s Bool -> Stmt s s  
assert e = \post -> let pre s = (e s) && (post s) in pre

{-
-- tougher. Needs loop invariant
while :: Expr s Bool -> Stmt s s -> Stmt s s
while e stmt = \post -> let pre s = if (e s) then ((while e stmt) (stmt post)) s else  in pre
-}

Finally here is a combination of the two above that uses the branching structure of the concrete execution to aid construction of the precondition. Although I haven’t expanded it out, we are using the full s t a b parametrization of lens in the sense that states go forward and predicates come back.


type Lens' a b = a -> (b, b -> a)
set l s a = let (_, f) = l s in f a


type Expr s a = s -> a
type Var s a = Lens' s a
type Pred a = a -> Bool
type Stmt s s' = s -> (s', Pred s' -> Pred s) -- eh. Screw the newtype

skip :: Stmt s s
skip = \x -> (x, id)


sequence :: Stmt s s' -> Stmt s' s'' -> Stmt s s''
sequence f g =   \s -> let (s', j) = f s in
                       let (s'', j') = g s' in
                           (s'', j . j')
assign :: Var s a -> Expr s a -> Stmt s s
assign v e = \s -> (set v s (e s), \p -> \s -> p (set v s (e s)))

--if then else
ite :: Expr s Bool -> Stmt s s' -> Stmt s s' -> Stmt s s'
ite e stmt1 stmt2 = \s -> 
                    if (e s) 
                    then let (s', wp) = stmt1 s in
                         (s', \post -> \s -> (e s) && (wp post s))
                    else let (s', wp) = stmt2 s in
                            (s', \post -> \s -> (not (e s)) && (wp post s))

assert :: Pred s -> Stmt s s
assert p = \s -> (s, \post -> let pre s = (post s) && (p s) in pre)

while :: Expr s Bool -> Stmt s s -> Stmt s s
while e stmt = \s -> if e s then let (s' , wp) = (while e stmt) s in
                                 (s', \post -> let pre s'' = (post s'') && (wp post s'') in pre)   
                            else (s, \p -> p)

{-

-- declare and forget can change the size and shape of the state space.
-- These are heterogenous state commpands
declare :: Iso (s,Int) s' -> Int -> Stmt s s'   
declare iso defalt = (\s -> to iso (s, defalt), \p -> \s -> p $ to iso (s, defalt)) 

forget :: Lens' s s' -> Stmt s s' -- forgets a chunk of state

declare_bracket :: Iso (s,Int) s' -> Int ->  Stmt s' s' -> Stmt s s
declare_bracket iso defalt stmt = (declare iso default) . stmt . (forget (_1 . iso))

Neat. Useful? Me dunno.

Relational Algebra with Fancy Types

Last time, I tried to give a primer of relations and relational algebra using the Haskell type type Rel a b = [(a,b)]. In this post we’re going to look at these ideas from a slightly different angle. Instead of encoding relations using value level sets, we’ll encode relations in the type system. The Algebra of Programming Agda repo and the papers quoted therein are very relevant, so if you’re comfortable wading into those waters, give them a look. You can find my repo for fiddling here

At this point, depending on what you’ve seen before, you’re either thinking “Yeah, sure. That’s a thing.” or you’re thinking “How and why the hell would you do such a ridiculous thing.”

Most of this post will be about how, so let’s address why first:

  1. Examining relations in this style illuminates some constructions that appear around the Haskell ecosystem, particularly some peculiar fellows in the profunctor package.
  2. More syntactic approach to relations allows discussion of larger/infinite domains. The finite enumerations of the previous post is nice for simplicity, but it seems you can’t get far that way.
  3. Mostly because we can – It’s a fun game. Maybe a useful one? TBD.

With that out of the way, let’s go on to how.

Translating Functions to Relation GADTs

We will be using some Haskell extensions in this post, at the very least GADTs and DataKinds. For an introduction to GADTs and DataKinds, check out this blog post. DataKinds is an extension that reflects every data constructor of data types to a type constructor. Because there are values True and False there are corresponding types created'True and 'False. GADTs is an extension of the type definition mechanism of standard Haskell. They allow you to declare refined types for the constructors of your data and they infer those refined type when you pattern match out of the data as well, such that the whole process is kind of information preserving.

We will use the GADT extension to define relational datatypes with the kind

a -> b -> *
. That way it has a slot a for the “input” and b for the “output” of the relation. What will goes in these type slots will be DataKind lifted types like 'True, not ordinary Haskell types like Int. This is a divergence from from the uses of similar kinds you see in Category, Profunctor, or Arrow. We’re doing a more typelevel flavored thing than you’ll see in those libraries. What we’re doing is clearly a close brother of the singleton approach to dependently typed programming in Haskell.

Some examples are in order for what I mean. Here are two simple boolean functions, not and and defined in ordinary Haskell functions, and their equivalent GADT relation data type.

not True = False
not False = True

data Not a b where
    NotTF :: Not 'True 'False
    NotFT :: Not 'False 'True

and True True = True
and False _ = False
and _ False = False

data And a b where
    AndTT :: And '( 'True, 'True) 'True
    AndFU :: And '( 'False, a) 'False
    AndUF :: And '( a, 'False) 'False

You can already start to see how mechanical the correspondence between the ordinary function definition and our new fancy relation type. A function is often defined via cases. Each case corresponds to a new constructor of the relation and each pattern that occurs in that case is the pattern that appears in the GADT. Multiple arguments to the relations are encoded by uncurrying everything by default.

Any function calls that occur on the right hand side of a function definition becomes fields in the constructor of our relation. This includes recursive calls and external function calls. Here are some examples with a Peano style natural number data type.

data Nat = S Nat | Z

plus Z x = x
plus (S x) y = S (plus x y)

data Plus a b where
    PZ :: Plus '( 'Z, a) a
    PS :: Plus '( a,b) c -> Plus '( 'S a, b) c 

We can also define things that aren’t functions. Relations are a larger class of things than functions are, which is part of their utility. Here is a “less than equal” relation LTE.

data LTE a b where
    LTERefl :: LTE n n
    LTESucc :: LTE n m -> LTE n ('S m)

You can show that elements are in a particular relation by finding a value of that relational type. Is ([4,7], 11) in the relation Plus? Yes, and I can show it with with the value PS (PS (PS (PS PZ))) :: Plus (4,7) 11 . This is very much the Curry-Howard correspondence. The type R a b corresponds to the proposition/question is (a,b) \in R .

The Fun Stuff : Relational Combinators

While you need to build some primitive relations using new data type definitions, others can be built using relational combinators. If you avoid defining too many primitive relations like the above and build them out of combinators, you expose a rich high level manipulation algebra. Otherwise you are stuck in the pattern matching dreck. We are traveling down the same road we did in the previous post, so look there for less confusing explanations of the relational underpinnings of these constructions, or better yet some of the references below.

Higher order relational operators take in a type parameters of kind

a -> b -> *
and produce new types of a similar kind. The types appearing in these combinators is the AST of our relational algebra language.

The first two combinators of interest is the composition operator and the identity relation. An element (a,c) is in R \cdot Q if there exists a b such that (a,b) \in R and (b,c) \in Q. The fairly direct translation of this into a type is

{- rcompose :: Rel b c -> Rel a b -> Rel a c  -}

data RCompose k1 k2 a c where
   RCompose :: k1 b c -> k2 a b -> RCompose k1 k2 a c

type k <<< k' = RCompose k k' 
type k >>> k' = RCompose k' k

The type of the composition is the same as that of Profunctor composition found in the profunctors package.

type RCompose = Procompose

Alongside a composition operator, it is a knee jerk to look for an identity relation and we do have one

data Id a b where
   Refl :: Id a a


-- monomorphic identity. Leave this out?
data IdBool a b where
  ReflTrue :: IdBool 'True 'True
  ReflFalse :: IdBool 'False 'False

This is also a familiar friend. The identity relation in this language is the Equality type.

-- identity function is the same as Equality
type Id a b = Id (a :~: b)

We can build an algebra for handling product and sum types by defining the appropriate relational combinators. These are very similar to the combinators in the Control.Arrow package.

-- Product types

data Fan k k' a b where
    Fan :: k a b -> k' a c -> Fan k k' a '(b,c)

type k &&& k' = Fan k k'

data Fst a b where
    Prj1 :: Fst '(a, b) a

data Snd a b where
    Prj2 :: Snd '(a, b) b

-- Sum type

data Split k k' a b where
    CaseLeft :: k a c -> Split k k' ('Left a) c
    CaseRight :: k' b c -> Split k k' ('Right b) c

type k ||| k' = Split k k'

data Inj1 a b where
    Inj1 :: Inj1 a ('Left a)
data Inj2 a b where
    Inj2 :: Inj2 a ('Right a)

-- some derived combinators
type Par f g = Fan (f <<< Fst) (g <<< Snd)
type Dup  = Fan Id Id
type Swap = Fan Snd Fst

The converse of relations is very interesting operation and is the point where relations really differ from functions. Inverting a function is tough. Conversing a relation always works. This data type has no analog in profunctor to my knowledge and probably shouldn't.

data RConverse k a b where
    RConverse :: k a b -> RConverse k b a
-- Shorter synonym
type RCon = RConverse

Relations do not have a notion of currying. The closest thing they have is

data Trans k a b where
    Trans :: k '(a,b) c -> Trans k a '(b,c)

Lattice Operators

For my purposes, lattices are descriptions of sets that trade away descriptive power for efficiency. So most operations you'd perform on sets have an analog in the lattice structure, but it isn't a perfect matching and you're forced into approximation. It is nice to have the way you perform these approximation be principled, so that you can know at the end of your analysis whether you've actually really shown anything or not about the actual sets in question.

? No. No... Yes? Oh. OH! IT IS!

The top relation holds all values. This is represented by making no conditions on the type parameters. They are completely phantom.

newtype Top a b = Top ()

Bottom is a relation with no inhabitants.

newtype Bottom a b = Bottom Void

The meet is basically the intersection of the relations, the join is basically the union.

newtype RMeet k k' a b = RMeet (k a b, k' a b)
type k /\ k' = RMeet k k'  

newtype RJoin k k' a b = RJoin (Either (k a b) (k' a b))
type k \/ k' = RJoin k k' 

A Lattice has an order on it. This order is given by relational inclusion. This is the same as the :-> combinator can be found in the profunctors package.

type (:->) p q = forall a b. p a b -> q a b
type RSub p q = p :-> q

Relational equality can be written as back and forth inclusion, a natural isomorphism between the relations. There is also an interesting indirect form.

data REq k k' = REq {to' :: k :-> k', from' :: k' :-> k }

Relational Division

If we consider the equation (r <<< p) :-> q with p and q given, in what sense is there a solution for r? By analogy, this looks rather like r*p = q, so we're asking a kind of division question. Well, unfortunately, this equation may not necessarily have a solution (neither do linear algebraic equations for that matter), but we can ask for the best under approximation instead. This is the operation of relational division. It also appears in the profunctor package as the right Kan Extension. You'll also find the universal property of the right division under the name curryRan and uncurryRan in that module.

newtype Ran p q a b = Ran { runRan :: forall x. p x a -> q x b }
type RDiv = Ran

One formulation of Galois connections can be found in the adjunctions file. Galois Connections are very slick, but I'm running out of steam, so let's leave that one for another day.

Properties and Proofs

We can prove many properties about these relational operations. Here a a random smattering that we showed using quickcheck last time.

prop_ridleft ::  (k <<< Id) :-> k
prop_ridleft (RCompose k IdRefl) = k

prop_ridright ::  (Id <<< k) :-> k
prop_ridright (RCompose IdRefl k) = k

prop_meet :: p /\ q :-> p
prop_meet (RMeet (p, q)) = p

prop_join :: p :-> p \/ q
prop_join p = RJoin (Left p)

meet_assoc :: RMeet k (RMeet k' k'') a b -> RMeet (RMeet k k') k'' a b
meet_assoc (RMeet (k, (RMeet (k',k'')))) = RMeet (RMeet (k,k'), k'')

prop_top :: k :-> Top
prop_top _ = top

prop_bottom :: Bottom :-> k
prop_bottom (Bottom x) = absurd x

bottom_compose :: REq (k <<< Bottom) Bottom
bottom_compose = REq (\(RCompose k (Bottom b)) -> absurd b) prop_bottom

data Iso a b = Iso {to :: a -> b, from :: b -> a}
type a <-> b = Iso a b

meet_universal :: (p ::-> RMeet k k') <-> (p ::-> k, p ::-> k')
meet_universal = Iso to from where
    to (RSub f) = (RSub $ \p -> case f p of RMeet (k,k') -> k  , RSub $ \p -> case f p of RMeet (k,k') -> k')
    from (RSub f,RSub g) = RSub $ \p -> RMeet (f p, g p) 

prop_con :: RCon (RCon k) :-> k
prop_con (RConverse (RConverse k)) = k

Odds and Ends

  • Recursion Schemes - Recursion schemes are a methodology to talk about recursion in a point free style and where the rubber meets the road in the algebra of programming. Here is an excellent series of articles about them. Here is a sample of how I think they go:
data MapMaybe k a b where
    MapJust :: k a b -> MapMaybe k ('Just a) ('Just b)
    MapNothing :: MapMaybe k 'Nothing 'Nothing

data Cata map k a b where
    Cata :: k fa a -> map (Cata map k) x fa  -> Cata map k ('Fix x) 
  • Higher Order Relations?
  • Examples of use. Check out the examples folder in the AoP Agda repo. These are probably translatable into Haskell.
  • Interfacing with Singletons. Singletonized functions are a specialized case or relations. Something like?
  •  newtype SFun a b = SFun (Sing a -> Sing b) 
  • A comment to help avoid confusion. What we've done here feels confusingly similar to profunctor, but it is in fact distinct I think. Profunctors are described as a categorical generalization of relations , but to be honest, I kind of don't get it. Despite many of our constructions appearing in the profunctor package, the profunctor typeclass itself appears to not play a role in our formulation. There just isn't a good way to dimap under our relations as written, unless you construct free profunctors. Converse at the least is a wrench in the works.
  • Star and graphs. Transition relations are a powerful methodology. A transition relation is in some respects the analog of a square matrix. We can iteratively compose it with itself.
-- Check out "term rewriting and all that"
-- This is also the reflection without remorse data type
-- TSequence http://okmij.org/ftp/Haskell/zseq.pdf
-- this is also a free instance of Category
data Star k a b where
    Done :: Star k a a
    Roll :: k b c -> Star k a b -> Star k a c

data KPlus k a b where
    PDone :: k a b -> KPlus k a b
    PRoll :: k b c -> KPlus k a b -> KPlus k a c

type SymClos k a b = RJoin k (RCon k) a b
type RefClos k a b = RJoin k Id a b
{- n-fold composition -}
-- similar to Fin.
-- This is also the Vec n is to list and this is to reflection without remorse. Kind of interesting
data NFold n k a b where
    One :: k a b -> NFold ('S n) k a b
    More :: k b c -> NFold n k a b -> NFold ('S n) k a b

References

CAV 2019 Notes: Probably Nothin Interestin’ for You. A bit of noodling with Liquid Haskell

I went to the opening workshops of CAV 2019 in New York this year (on my own dime mind you!) after getting back from joining Ben on the long trail for a bit. The tutorials on Rosette and Liquid Haskell were really fun. Very interesting technology. Also got some good ramen and mochi pudding, so it’s all good. Big Gay Ice Cream was dece.

Day 1 workshops

Calin Belta http://sites.bu.edu/hyness/calin/.Has a new book. Control of Temporal logic systems. Automata. Optimized. Partition space into abstraction. Bisimulation https://www.springer.com/gp/book/9783319507620

Control Lyapunov Function (CLF) – guarantees you are going where you want to go

Control Barrier Function – Somehow controls regions you don’t want to go to.

Lyapunov function based trajectory optimization. You somehow have (Ames 2014) http://ames.gatech.edu/CLF_QP_ACC_final.pdf Is this it?

Differential flatness , input output linearization

Sadradiini worked there.

Temproal logic with

Rise of Temporal Logic

Linear Temporal Logic vs CTL

Fixpoint logic,

Buchi automata – visit accepting state infinite times

equivalency to first order logic

monadic logic, propositions only take 1 agrument. Decidable. Lowenheim. Quantifier elimination. Bounded Mondel property

Languages: ForSpec, SVA, LDL, PSL, Sugar

Monadic second order logic (MSO).

method of tableau

Sadraddini

Polytopic regions. Can push forward the dynmaics around a trajectory and the polytope that you lie in. RRT/LQR polytopic tree. pick random poitn. Run.

Evauating branching heuristics

branch and prune icp. dreal.

branch and prune. Take set. Propagate constraints until none fire.

branching heuristics on variables

largest first, smearing, lookahead. Try different options see who has the most pruning. Non clear that helped that muhc

QF_NRA. dreal benchmarks. flyspeck, control, robotics, SMT-lib

http://capd.sourceforge.net/capdDynSys/docs/html/index.html

Rosette

commands: saolver adied programming

verify – find an input on which the assertions fail. exists x. not safe

debug – Minimal unsat core if you give an unsat query. x=42/\ safe(s,P(x))$ we know thia is unsat because of previous step

solve – exists v si.t safe(v)

synthesis – exists e forall x safe(x,P(x))

define-symbolic, assert, verify, debug, solve, sythesis

Rosette. Alloy is also connected to her. Z Method. Is related to relational logic?

https://homes.cs.washington.edu/~emina/media/cav19-tutorial/index.html

http://emina.github.io/rosette/

Building solver aided programming tool.

symbolic compiler. reduce program all possible paths to a constraint

Cling – symbolic execution engine for llvm

implement intepreter in rosette

Symbolic virtual machine

layering of languages. DSL. library (shallow) embedding. interpreter (deep) embedding.

deep embedding for sythesis.

I can extract coq to rosette?

how does it work?

reverse and filter keeping only positive queries.

symbolic execution vs bounded model checking

symbolic checks every possible branch of the program. Cost is expoentntial

CBMC.

type driven state merging. Merge instances of primitiv types. (like BMC), value types structurally ()

instance Merge Int, Bool, Real — collect up SMT context

vs. Traversable f => Merge (f c) – do using Traversable

symbolic union a set of guarded values with diskoint guard.

merging union. at most one of any shape. bounded by number of possible shapes.

puts some branching in rosette and some branch (on primitives) in SMT.

symbolic propfiling. Repair the encdoing.

tools people have built.

veify radiation

strategy generation. That’s interesting. builds good rewrite rules.

serval.

certikso komodo keystone. fintie programss

IS rosette going to be useful for my work? coooooould be

Liquid Haskell

https://ranjitjhala.github.io/

{-# LANGUAGE GADTs, DataKinds, PolyKinds #-}
{-@ LIQUID "--reflection" @-}
{-@ LIQUID "--short-names" @-}
{-@ LIQUID "--ple" @-}

{-@ type TRUE = {v: Bool | v = True }@-}
{-@  type NonEmpty a = {v : [a] | len v > 0}  @-}
{-@ head :: NonEmpty a -> a @-}
head (x : _) = x

{-@ measure f :: Int -> Int @-}
f x = 2 * x

{-@ true :: TRUE @-}
true = True


-- impl x y = x ==> y
-- {-@ congruence :: Int -> Int -> TRUE @-}
-- congruence x y = (x == y) ==> (f x == f y)
-- {-@ (==>) :: {x : Bool |} -> {y : Bool |} -> {z : Bool | z = x ==> y} @-}
-- x ==> y = (not x) || y

-- aws automated reaosning group
{-
nullaway uber. 
sadowski google static anaylis
infer faecobook

give programmer early 


refinement types

why and how to use
how to implement refinement types

mostly like floyd hoare logic

types + predicates = refinement type

t := {x:b | p} erefinemed b
      x : t -> t -- refined function type
linear arithemtic

congruence axioms


emit a bunch of verification conditions (VC)
p1 => p2 => p3 ... 

SMT can tell if VC is alwasy true


-}

{-@ type Zero = {v : Int | v == 0} @-}
{-@ zero :: Zero @-}
zero = (0 :: Int) -- why?


{-@ type NAT = {v: Int | v >= 0} @-}
{-@ nats :: [Nat]  @-}
nats = [0,1,1,1] :: [Int]

{-

subtype

in an environemnt Gamma, t1 is a subtype of t2
x are vairables are in scope.
/\x |- {v | q} <= {v | r}

True  => 


-}


{-@ plus :: x : Int -> y : Int -> {v : Int | v = x + y} @-}
plus :: Int -> Int -> Int
plus x y = x + y

-- measure. uninterpeteyd function called measure.
-- {-@ measure vlen :: Vector a -> Int @-}

-- {-@ at :: vec :Vector a -> {i : Nat | i < vlen vec} ->  @-}

-- {-@  size :: vec : Vector a -> {n : Nat | n == vlen vec} @-}

{-
horn contrints infer 

Collections and Higher and order fucntions

reduce :: (a -> b -> a) -> a -> [b] -> a
type a is an invaraint that holds on initial acc and indictively on f
Huh. That is true.


-- Huh. I could prove sum formulas this way.

sum'' = vec = let is = range 0 (size vec)
              add = \tot i -> ot + at vec i


              properties of data structures


size of of a list
allow uninterpetyed functions inside refinements

{ measre length}

LISTNE a = {v : [a] | 0 < length v}


measure yields refined constructions

[] :: {v : [a] | legnth v = 0}


              --}


              {-
              
              Q: measure?
              Q : environment?
              Q : where is standard libraru?
              
              no foralls anywhere? All in decidable fragment


              p : ([a],[a]) | (fst p) + (snd p) == lenght xs
              measures fst and snd


              interpeter

              impossible :: {v : String | false} -> a

imperative language

interpolation


/inlcude folder has prelude

basic values
{v : Int | lo <= v && v < hi }


invaraint properties of structures

encoe invraints in constructor

data OrdPair = OP {opX :: Int, opY :: Int}

{-@ data OrdPair = OP {opX :: Int, opY :: {v : Int | opX < v}@-}


class Liquid Int
class Liquid Bool
class Liquid ... Real?

{-@   @-}

Liquid Relations?
{  }

data OList 
{-@data OList a LB = Empt | :< {oHd :: {v : a | LB = v}   oTl :: OList a oHd }   @-}


{-@   {oHd :: a, oTl :: OList {v : a | oHd < v}   }@-}
              
GADTs?
-}
data MyList a where
    Nil :: MyList a
    {-@ Cons :: v : a -> MyList {x : a | x < v} ->  MyList a @-}
    Cons :: a -> MyList a -> MyList a

test :: MyList Int
test = Cons 2 (Cons 1 Nil)

{-

abstracting the invaraint from the data structure
parametrzie by relations

data [a]<rel :: a -> a -> Bool> where
    = []
    | (:) {hd :: }

    rel != is unique list

{\x y -> x >= y}
type level lambdas!? .... uh.... maybe.

reflecting singletons into liquid?

termination metrics

/ [length xs + len ys] -- merge sort

{-@ Half a s = }@-}

Oncey ou have temrination proofs you have proofs of correctness

Propositions as Types

Plus commutes is trivial
{n = n + n}

-}


{-@ easyProof :: {True} @-}
easyProof = ()

-- hot damn. I mean this is in it's legerdomain. But prettttty sweet.
{-@ commute :: x : Int -> y : Int -> {x + y = y + x} @-}
commute :: Int -> Int -> ()
commute x y = ()

{-@ reflect mysum @-}
{-@ mysum :: Nat -> Nat @-}
mysum :: Int -> Int
mysum 0 = 0 -- if n <= 0 then 0 else  2 * n + (mysum (n - 1))
mysum n = 2 * n + (mysum (n - 1))

-- what is going on here? why do I need _?
{-@ mysumpf :: _ -> {mysum 0 = 0 } @-}
-- mysumpf :: Proof
mysumpf _ = let x = mysum 0 in x

{-@ mysumpf' :: {mysum 3 = 12 } @-}
-- mysumpf :: Proof
mysumpf' = ()

{-@ reflect fastsum @-}
{-@ fastsum :: Nat -> Nat @-}
fastsum :: Int -> Int
fastsum n = n * (n + 1)

type Proof = ()
{-
{-@ pfsum :: x : Nat -> {fastsum x = mysum x} @-}
pfsum :: Int -> Proof
pfsum 0 = () -- let _ = fastsum 0 in let _ = mysum 0 in ()
pfsum n = pfsum (n-1)  
-}

{-@ pfsum :: x : Nat -> {fastsum x = mysum x} @-}
pfsum :: Int -> Proof
pfsum 0 = () -- let _ = fastsum 0 in let _ = mysum 0 in ()
pfsum n = pfsum (n-1)  

{-

reflection

reflect takes the prcondition of sum and dumps it as the poscondition

sum3 _ = let s0 =sum 0
             s1  = sum 1
             s2  = sum 3 -- all are going to be in scope.
z3 will connect the dots.



using proof combinatos from Proof Combinators
long chains of claculations

reflection of singletons

data SS s where
    {-@ SZero :: {v : Int | v = 0} -> SS 'Zero @-} 
    SZero :: Int -> SS 'Zero 
    {-@ SZero :: {v : Int | v = 0} -> SS 'S a @-} 
    SZero :: Int -> SS 'Zero 





    proof by induction
    sum n = n * (n + 1)/2

    2 * sum n  = n * (n + 1)


point free liquid types

(.) :: (a -> b) -> (a -> )
? Can I abstract over predicates like this?
({v:a | p} -> {s:}) -> 



Vectors

cauchy schwartz



-}

data V2 a = V2 a a 

{-@ reflect dot @-} 
dot (V2 x y) (V2 x' y') = x * x' + y * y'

{-@ reflect vplus @-}
vplus (V2 x y) (V2 x' y') = V2 (x + x') (y + y')

{-@ reflect smul @-}
smul s (V2 x' y') = V2 (s * x') (s * y')
{-
{-@ cauchy :: x : V2 Int -> y : V2 Int -> {(dot x y) * (dot x y) <= (dot x x) * (dot y y) } @-}
cauchy :: V2 Int -> V2 Int -> Proof
cauchy x y = let q = dotpos (vplus x y) in let r = dotpos (vplus x (smul (-1 :: Int) y)) in (\_ _ -> ()) q r
-}

-- {-@ square :: Int -> Nat @-} -- basiclly the same thing
{-@ reflect square @-}
square :: Int -> Int
square x = x * x



{-@ sqpos :: x: Int -> {square x >= 0} @-}
sqpos :: Int -> ()
sqpos x = ()

{-@ dotpos :: x: V2 Int -> {dot x x >= 0} @-}
dotpos :: V2 Int -> ()
dotpos x = ()

{-@ dotsym :: x: V2 Int -> y : V2 Int -> {dot x y = dot y x} @-}
dotsym :: V2 Int -> V2 Int ->  ()
dotsym x y = ()

{-@ vpluscomm :: x: V2 Int -> y : V2 Int -> {vplus x y = vplus y x} @-}
vpluscomm :: V2 Int -> V2 Int ->  ()
vpluscomm x y = ()

{-@ dotlin :: x: V2 Int -> y : V2 Int -> z : V2 Int -> {dot (vplus x y) z = dot x z + dot y z} @-}
dotlin :: V2 Int -> V2 Int ->  V2 Int ->  ()
dotlin x y z = ()





{-

What else is interesting to prove?

verify stuff about ODEs?

fold [1 .. t] where t = 10

could give little spiel about how dynamical systems are like imperative programming



get some rationals.

profunctor
p a b
a -> b are refined functions


I should learn how to abstract over typeclasses. Verified typeclasses?

SMT has built in rationals prob?
-}

data Rat = Rat Int Int 

{-@ reflect rplus @-}
rplus :: Rat -> Rat -> Rat
rplus (Rat x y) (Rat x' y') = Rat (x*y' + x'*y) (y * y')



{-@ reflect rmul @-}
rmul :: Rat -> Rat -> Rat
rmul (Rat x y) (Rat x' y') = Rat (x*x') (y * y')





data Nat' = S Nat' | Z

{-@ measure nlen @-}
{-@ nlen :: Nat' -> Nat @-}
nlen :: Nat' -> Int
nlen Z = 0
nlen (S x) = 1 + (nlen x) 
{-
-- failing?
-- crash: SMTLIB2 respSat = Error "line 31 column 169: unknown sort 'Main.SNat'"
data SNat a where
    SZ :: SNat 'Z
    SS :: SNat x -> SNat ('S x)
-}

{-@ reflect conv @-}
{-@ conv :: x : Nat -> {v : Nat' | nlen v = x} @-}
conv :: Int -> Nat'
conv 0 = Z
conv x = S (conv (x-1))

-- It's an isomorphism

{-@ pfconv :: x : Nat -> {nlen (conv x) = x} @-}
pfconv :: Int -> Proof
pfconv 0 = ()
pfconv x = pfconv (x - 1)

{-@ pfconv' :: x : Nat' -> {conv (nlen x) = x} @-}
pfconv' :: Nat' -> Proof
pfconv' Z = ()
pfconv' (S x) = pfconv' x

{-@ reflect plus' @-}
plus' :: Nat' -> Nat' -> Nat'
plus' Z x = x
plus' (S x) y = S (plus' x y)

{-@ plusz' :: x : Nat' -> {plus' x Z = plus' Z x}  @-}
plusz' :: Nat' -> Proof
plusz' Z = ()
plusz' (S x) = plusz' x


{-@ pluscomm' :: x : Nat' -> y : Nat' -> {plus' x y = plus' y x} / [nlen x, nlen y] @-}
pluscomm' :: Nat' -> Nat' -> Proof
pluscomm' Z y = plusz' y
pluscomm' (S x) (S y) = const (pluscomm' (S x) y) $ const (pluscomm' x (S y)) $  pluscomm' x y
    -- const () $ const (plus' (S x) (S y)) $ const (plus' x (S y)) (plus' x y)
    
    -- const (pluscomm' (S x) y) $ const (pluscomm' x (S y)) $  pluscomm' x y
 -- flip const is proof combinator .==   
    {-let q = pluscomm' x (S y) in
                        let w = pluscomm' (S x) y in 
                        let r = pluscomm' x y in (\b n m -> ()) q w r -- ? Was this necessary? -}
pluscomm' x Z = plusz' x




-- {-@ data Iso = @-}
data Iso a b = Iso { to :: a -> b, from :: b -> a, p1 :: Proof, p2 :: Proof}


{-

We also have type level lambdas.
refinement polymorphism

LH is somewhat like singletons in the sense there is a manual reflection step.
In singletons the manual reflection is in the Sing type
in LH it is kind of all over the place. (+) has a type. Where is it defined?


How does it know that the Haskell function + is the same as the SMT solver function?

Coq and Agda and Idris type checking is powered quite a bit by an internal unification engine
explicit annotation may lessen the burden somewhat

SMT solvers as a unification engine
structure unification vs uninterpeted functions.
f a ~ Int is not a valid Haskell constraint. Maybe with the unmatchable arrow it is?
In a funny sense, there is a difference between  Just and (+ 1). 
One being a constructor means we can match out of it
Just :: a ->> b
(+ 1) :: Int -> Int

-}

-- test' :: (f a ~ Int) => ()
-- test' = ()

Liquid Haskell – What is?

another thing we could do is galois connections between refinements. Pos, Zero, Neg <-> Int

Liquid Haskell uses SMT solvers to resolve it’s type checking requirements.

Agda et al also work very much via unification. Unification is a broad term but it’s true.

It also has a horn clause solver for inference. Every language needs some kind of inference or you’d go insane. Also it is piggybacking on haskell

It’s not as magical as I thought? Like seeing the magicians trick. It really does understand haskell code. Like it isn’t interpretting it. When it knows facts about how (+) works, that is because the refined type was put in by hand in the prelude connecting it to SMT facts. What is imported by liquid haskell?

The typing environment is clutch. You need to realize what variables are in scope and what their types are, because that is all the SMT can use to push through type checking requirements.

Installing the stack build worked for me. It takes a while . I couldn’t get cabal install to work, because I am not l33t.

Uninterpeted functions. Unmatchability?

It wouldn’t be haskell without a bunch of compiler directives. It is somewhat difficult to find in a single cohesive place what all the syntax, and directives are from liquid haskell. Poking around it best.

  • ple
  • reflection
  • no-termination
  • higherorder – what is this?

https://github.com/ucsd-progsys/230-wi19-web course notes

https://github.com/ucsd-progsys/liquid-sf some of software foundations

https://nikivazou.github.io/publications.html niki vazou’s pubs. Check out refinement reflection

https://nikivazou.github.io/static/Haskell17/law-abiding-instances.pdf draft work? Shows stuff about typeclasses. This is a haskell 2017 paper though

https://arxiv.org/pdf/1701.03320 intro to liquid haskell. Interesting to a see a different author’s take

http://goto.ucsd.edu/~nvazou/presentations/ presentations. They are fairly similar to one another.

Liquid haskell gives us the ability to put types on stuff that wasn’t possible before.

Linearity :: f :: {a -> b | f (s ^* a) == s ^* (f a) }

Pullback. {(a,b) | f a == g b}

Equalizer

Many things in categoruy theory rely on the exists unique. Do we have functiona extensionality in Liquid haskell?

product : {(a,b) | f q = x, g q = y, => }

Pushing the boundaries on what liquid haskell can do sounds fun.

Equalizer. The eqaulizer seems prominents in sheaves. Pre-sheaves are basically functors. Sheaves require extra conditions. Restriction maps have to work? Open covers seem important

type Equalizer f g a b = {(e :: a , eq :: a -> b) | f (eq e) = g (eq e) }

I think both the type a and eq are special. e is like an explcit parametrization.

type Eq f g a = {e :: a | f e = g e} I think this is more in the spirit. Use f and g both as measures.

presheaf is functor. But then sheaf is functor that

(a, Eq (F a) (G a)). typelevel equalizer? All types a that F and G agree on.

https://ncatlab.org/nlab/show/equalizer

https://blog.functorial.com/posts/2012-02-19-What-If-Haskell-Had-Equalizers.html

Records are sheaves – Jon Sterling. Records have subtyping. This gives you a toplogy feeling thing.

https://www.slideshare.net/jonsterling/galois-tech-talk-vinyl-records-in-haskell-and-type-theory

What about purescript records?

{foo | a} {bar | a} -> intersection = {foo bar | b} can inhabit either

union is
or do you want closed records? union is union of fields. intersection is intersection of fields.

In this case a cover would be a set of records with possibly overlapping fields whose combined labels cover the whle space we want to talk about. consistency condition of sheaf/equalizer is that overlapping records fields have to match. I guess {  q.foo = r.foo } ?There is a way to combine all the stuff up. This is exactly what Ghrist was getting at with tables. Tables with shared columns.

data R1 = R1 {foo :: Int, bar :: Int}

{ (r1 :: R1, r2 :: R2) | (foo r1) = (foo r2) } — we manitain duplicates across records.

{. }

if you have a “cover” {foo bar |} {bar fred} {gary larry} whose in

https://www.sciencedirect.com/science/article/pii/S1571066108005264

Sheaves. As a model of concurrency? Gaguen paper.

sheaves as constraint satisfcation? sheafifcation. Constraint solving as a way of fusing the local constraints to be globally consistent.

sheaves as records

sheaves as data fusion

http://www.cs.bham.ac.uk/~mhe/papers/barbados.pdf

Escardo. Compact data types are those finitely searchable

Continuous funcitons are ~computable? Productive?

http://www.paultaylor.eu/

http://www.paultaylor.eu/ASD/foufct/

http://www.paultaylor.eu/~pt/prafm/

typed recursion theory toplogy

typed computatabiltity theory

Topological notions in computation. Dictionary of terms realted decidable, searchable, semi decidablee

cs.ioc.ee/ewscs/2012/escardo/slides.pdf

https://en.wikipedia.org/wiki/Computable_topology

Through NDArray overloading, a significant fragment of numpy code is probably verifiable.

Start with functional arithmetic programs.

Need to inspect function annotations to know how to build input type.

@verify() tag

Use (Writer a) style monad.

If statements are branching. We are again approaching inspecting functions via probing. But what if we lazily probe. At every __bool__ point, we run a z3 program to determine if there is an avaiable bool left possible (we don’t need to inspect dead code regions. Also would be cool to mention it is a dead region). Curious. We’re reflecting via Z3.

Loops present a problem. Fixed loops are fine. but what about loops that depend on the execution? for i in range(n). I guess again we can hack it…? Maybe. range only takes an integer. we don’t have overload access.

Maybe we need to go into a full eval loop. utterly deconstructing the function and evaluating it statelemnt by statement.

(compare :: a -> a -> Comparison). We could select a choice based on if there is a new one avaialable. Requires access to external store. We lose the thread. How can we know a choice was made? How can we know what the choice was? Did it ask var1 or var2? We can probably do it in python via access to a global store. But in haskell?

while loops take invariant annotations.

It would be cool to have a program that takes

pre conditions. Post ocnditions, but then also a Parameter keyword to declare const variables as deriveable. exists parameter. forall x precondition x => post condition.

Parameter could be of a type to take a dsl of reasonable computations. Perhaps with complexity predicates. and then interpretting the parameter defines the computation.

Or simpler case is parameter is an integer. a magic number.



@pre(lambda x: None)
@post(lambda r: r >= 0)
def square(x):
	return x**2

@verify(pre, post) # Easier. because now we can also verify the individual function. Call Z3 at function definition time.
def pre(f,cond):
   if(VERIFCAIOTN_ON)
   return fnew
   def fnew(x):
        if(VERIFICATION_ON):
	   		if(x == VerificationEnv):
	   			newenv = x.copy
	            new.add_pre(cond(x.var))
	            newVar = Z3.variable()
	            newenv.add(newVar == f(x.var))



	        else:
	        	return f(x)



def post(f, cond):
   def fnew(x):
    if x == VerifcationEnv:
        Z3.findmodel(not cond(x.var), x.env)
        if can't find one:
           we're good.
           x.env.append(cond(x.var))
           return x
        else:
           assert(False, model, function name, function postcondition code.)


#overloading assignment. isn't a problem.
class VerifyArray():
   #numpy Z3 shim.

#termination requires measure decreasing at every recursive call.
#arbaitrary loops? How to deal with those?
#maybe a hierarchy of envs to make it easier on z3. Like it doesn't need to give the whole program history if the local use is good enough.



class VerificationEnv():
    self.var = []
    self.pre = []
    self.post = []

Proving some Inductive Facts about Lists using Z3 python

Z3 is an SMT solver which has a good rep. Here are some excellent tutorials.

https://rise4fun.com/z3/tutorial

https://theory.stanford.edu/~nikolaj/programmingz3.html

http://homepage.cs.uiowa.edu/~ajreynol/VTSA2017/

SMT stands for satisfiability modulo theories. The exact nature of power of these kinds of solvers has been and is still hazy to me. I have known for a long time that they can slam sudoku or picross or other puzzles, but what about more infinite or logic looking things? I think I may always be hazy, as one can learn and discover more and more encoding tricks to get problems and proofs that you previously thought weren’t solvable into the system. It’s very similar to learning how to encode to linear programming solvers in that respect.

SMT solvers combine a general smart search procedure with a ton of specialized solvers for particular domains, like linear algebra, polynomials, linear inequalities and more.

The search procedure goes by the name DPLL(T). It is an adjustment of the procedure of SAT solvers, which are very powerful and fast. SAT solvers find an assignment of boolean variables in order to make a complicated boolean expression true, or to eventually find that it can never be made true. SAT solvers work on the principle of guessing and deducing. If a OR b needs to be true and we already know a is false, we can deduce b must be true. When the deduction process stalls, we just guess and then backtrack if it doesn’t work out. This is the same process you use manually in solving Sudoku.

The modern era of SAT solvers came when a couple new tricks vastly extended their power. In particular Conflict Driven Clause Learning (CDCL), where when the solver finds itself in a dead end, it figures out the choices it made that put it in the dead end and adds a clause to its formula so that it will never make those choices again.

https://sahandsaba.com/understanding-sat-by-implementing-a-simple-sat-solver-in-python.html

SMT works by now having the boolean variables of the SAT solver contain inner structure, like the boolean p actually represents the fact x + y < 5. During the search process it can take the pile of booleans that have been set to true and ask a solver (maybe a linear programming solver in this case) whether those facts can all be made true in the underlying theory. This is an extra spice on top of the SAT search.

Something that wasn’t apparent to me at first is how important the theory of uninterpreted formulas is to SMT solvers. It really is their bread and butter. This theory is basically solved by unification, which is the fairly intuitive process of finding assignments to variables to make a set of equations true. If I ask how to make fred(x,4) = fred(7,y), obviously the answer is y=4, x=7. That is unification. Unification is a completely syntax driven way to deduce facts. It starts to give you something quite similar to first order logic.

https://eli.thegreenplace.net/2018/unification/

https://cs.mtsu.edu/~rbutler/courses/sdd/automated_deduction/unification.pdf

I was also under the impression that quantifiers \forall, \exists were available but heavily frowned upon. I don’t think this is quite true. I think they are sort of a part of the entire point of the SMT solver now, although yes, they are rather flaky. There are a number of methods to deal with the quantifier, but one common one is to look for a pattern or parts of the subformula, and instantiate a new set of free variables for all of the quantified ones and add the theorem every time the patterns match. This is called E-matching.

Here are a couple tutorials on proving inductive facts in SMT solvers. You kind of have to hold their hand a bit.

http://lara.epfl.ch/~reynolds/vmcai15.pdf

http://homepage.divms.uiowa.edu/~ajreynol/pres-ind15.pdf

SMT solvers queries usually have the flavor of finding something, in this case a counterexample. The idea is that you try to ask for the first counterexample where induction failed. Assuming that proposition P was true for (n-1), find n such that P is not true. If you can’t find it, then the induction has gone through.

And here is some code where I believe I’m showing that some simple list functions like reverse, append, and length have some simple properties like \forall t. rev (rev(t)) == t .

from z3 import *

# Very useful reference
# https://theory.stanford.edu/~nikolaj/programmingz3.html

f = Function('f', IntSort(), IntSort())
s = Solver()
#s.add(f(True) == False, f(False) == True)
x = Int('x')
s.add(ForAll([x], f(x) >= x)) #> and < do not seem to be returning
print(s.sexpr())
s.check()
print(s.model())

# Rolling my own list data type
# Z3 has a built in which will probably be better?
s = Solver()

L = Datatype("MyList")
L.declare("Nil")
L.declare("Cons", ("hd", IntSort()), ("tail", L))
L = L.create()

t = Const('t', L)
u = Const('u', L)

y = Int('y')


rev = Function('reverse', L, L)
app = Function('append', L, L, L)
leng = Function('length', L, IntSort())

#defining my functions. Micro Haskell, BABY
#length
s.add( leng(L.Nil) == 0 )
s.add(  ForAll([u,y],  leng(L.Cons(y,u)) == 1 + leng(u))) #  patterns = [leng(L.Cons(y,u))] 

#append
s.add(  ForAll([u], app(L.Nil, u) == u)) 
s.add(  ForAll([t, u, y] , app(L.Cons(y,t), u)  ==  L.Cons(y, app(t, u ))))

#reverse
s.add( rev(L.Nil) == L.Nil)
s.add(  ForAll([y,t],  rev(L.Cons(y,t)) == app(rev(t), L.Cons(y, L.Nil))))


s.push()
print("proving leng(t) >= 0")
#s.add( Or(And(t == L.Cons(y,u),  leng(u) >= 0 ), t == L.Nil))
s.add(  Not(leng(t) >= 0 ))
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),  leng(L.tail(t)) >= 0 ))
print(s.check())

s.pop()
s.push()
#s.add( leng(app(L.Nil, u)) == leng(u) )

print("prove length is preserved under app.")
s.add( leng(app(t,u)) != leng(t) + leng(u))
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   leng(app(L.tail(t),u)) == leng(L.tail(t)) + leng(u)  ))
print(s.check())

s.pop()
s.push()

print("reverse preserves length")
#Lemma Could place in clause with the above proof of this theorem
s.add( ForAll([t,u], leng(app(t,u)) == leng(t) + leng(u)   )) #subgoal needed
s.add( leng(rev(t)) != leng(t))
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   leng(rev(L.tail(t))) == leng(L.tail(t)) ))


s.pop()
s.push()

print("reverse reverse = id")
s.add( ForAll( [t,u], rev(app(t,u)) ==  app(rev(u), rev(t)) ) ) #subgoal needed
s.add( rev(rev(t)) != t )
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   rev(rev(L.tail(t))) == L.tail(t) ))
print(s.check())


#junk

#s.add(t != L.Nil ) 
#s.add( ForAll([t], rev(L.Cons(y,t)) ==   ) , rev(L.Nil) == L.Nil)
#s.add( leng(L.Cons()) == 0 )
#s.add(  ForAll(y,  leng(L.Cons(y,L.Nil)) == 1 + leng(L.Nil)))
#s.add(  Not( leng(app(t,u))  == leng(t) + leng(u) )) 

# prove length of app + nil is same
#s.add( leng(app(t,L.Nil)) != leng(t))
#s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   leng(app(L.tail(t),L.Nil)) == leng(L.tail(t))))

#s.add( app(L.Nil,L.Nil) == L.Nil)
#s.add( app(t,L.Nil) != t)
#s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   app(L.tail(t),L.Nil) == L.tail(t)))


#s.add(Or(t == L.Nil , And( t == L.Cons(y, u),   app(u,L.Nil) == u)  ))
#s.add( Implies(t == L.Nil, leng(app(L.Nil, u)) == leng(u) ))
# all of these don't work
#s.add( Implies(u == L.tail(t),  leng(u) >= 0 ))
#s.add( ForAll(u, Implies(t == L.Cons(y,u),  leng(u) >= 0 )))
#print(s.sexpr())
#print(s.check())
#print(s.model())

#print(tactics())
'''
def induction(freevar, f, construtors?):
    x = Int('x')
    return And(Not(f(x)), 
'''

'''
# controller synthesis
pos = Array(Real, 10) 
for t in range(10):
    s.add(pos[t] == pos[t] +  v * dt)
    s.add(v[t] == v[t] + f(pos[t]) * dt)
    s.add(f(pos[t] <= 1)) # constraints on force
s.add(ForAll([init], Implies(And(init <= 1, pos[0] == init), pos[9] <= 0.1) ))
 
s.add(Forall[x], f(x) == a * x) # parametrizing. 
'''
#s.set("produce-proofs", True
#s.set("mbqi", False)
#s.set(mbqi=False)
#print(s.proof())