Casadi – Pretty Damn Slick

Casadi is something I’ve been aware of and not really explored much. It is a C++ / python / matlab library for modelling optimization problems for optimal control with bindings to IPOpt and other solvers. It can produce C code and has differentiation stuff. See below for some examples after I ramble.

I’ve enjoyed cvxpy, but cvxpy is designed specifically for convex problems, of which many control problems are not.

Casadi gives you a nonlinear modelling language and easy access to IPOpt, an interior point solver that works pretty good (along with some other solvers, many of which are proprietary however).

While the documentation visually looks very slick I actually found it rather confusing in contents at first. I’m not sure why. Something is off.

You should download the “example pack” folder. Why they don’t have these in html on the webpage is insane to me. https://github.com/casadi/casadi/releases/download/3.4.4/casadi-example_pack-v3.4.4.zip

It also has a bunch of helper classes for DAE building and other things. They honestly really put me off. The documentation is confusing enough that I am not convinced they give you much.

The integrator classes give you access to external smart ode solvers from the Sundials suite. They give you good methods for difficult odes and dae (differential algebraic equations, which are ODEs with weird constraints like x^1 + y^1 == 1) Not clear to me if you can plug those in to an optimization, other than by a shooting method.

Casadi can also output C which is pretty cool.

I kind of wondered about Casadi vs Sympy. Sympy has lot’s of general purpose symbolic abilities. Symbolic solving, polynomial smarts, even some differential equation understanding. There might be big dividends to using it. But it is a little harder to get going. I feel like there is an empty space for a mathemtical modelling language that uses sympy as it’s underlying representation. I guess monkey patching sympy expressions into casadi expression might not be so hard. Sympy can also output fast C code. Sympy doesn’t really have any support for sparseness that I know of.

As a side note, It can be useful to put these other languages into numpy if you need extended reshaping abilities. The other languages often stop at matrices, which is odd to me.

Hmm. Casadi actually does have access to mixed integer programs via bonmin (and commercial solvers). That’s interesting. Check out lotka volterra minlp example

https://groups.google.com/forum/#!topic/casadi-users/8xCHmP7UmpI

The optim interface makes some of this look better. optim.minimize and subject_to. Yeah, this is more similar to the interfaces I’m used to. It avoids the manual unpacking of the solution and the funky feel of making everything into implicit == 0 expressions.

https://web.casadi.org/blog/opti/

Here is a simple harmonic oscillator example using the more raw casadi interface. x is positive, v is velocity, u is a control force. I’m using a very basic leap frog integration. You tend to have to stack things into a single vector with vertcat when building the final problem.

from casadi import *
import matplotlib.pyplot as plt

g = 9.8
N = 100

x = SX.sym('x',N)
v = SX.sym('v', N)
u = SX.sym('u', N-1)
#theta = SX('theta', N)
#thdot = SX('thetadot', N)

dt = 0.1
constraints = [x[0]-1, v[0]] # expressions that must be zero
for i in range(N-1):
    constraints += [x[i+1] - (x[i] + dt * v[i]) ]
    constraints += [v[i+1] - (v[i] - dt * x[i+1] + u[i] * dt)]

cost = sum([x[i]*x[i] for i in range(N)]) + sum([u[i]*u[i] for i in range(N-1)])

nlp = {'x':vertcat(x,v,u), 'f':cost, 'g':vertcat(*constraints)}
S = nlpsol('S', 'ipopt', nlp)
r = S(lbg=0, ubg=0) # can also give initial solutiuon hint, some other things
x_opt = r['x']
x = x_opt[:N]
v = x_opt[N:2*N]
u = x_opt[2*N:]
#u_opt = r['u']
print('x_opt: ', x_opt)
print(S)
plt.plot(x)
plt.plot(u)
plt.plot(v)
plt.show()

Let’s use the opti interface, which is pretty slick. Here is a basic cartpole https://web.casadi.org/blog/opti/

from casadi import *
import matplotlib.pyplot as plt

g = 9.8
N = 100
# https://web.casadi.org/blog/opti/



opti = casadi.Opti()

x = opti.variable(N)
v = opti.variable(N)
theta = opti.variable(N)
dtheta = opti.variable(N)
u = opti.variable(N-1)


opti.subject_to( u <= 1) 
opti.subject_to( -1 <= u) 
opti.subject_to( x <= 2) 
opti.subject_to( -2 <= x) 
opti.subject_to(x[0] == 0)
opti.subject_to(v[0] == 0)
opti.subject_to(theta[0] == 0)
opti.subject_to(dtheta[0] == 0)

dt = 0.05
for i in range(N-1):
    opti.subject_to( x[i+1] == x[i] + dt * (v[i]))
    opti.subject_to( v[i+1] == v[i] + dt * (x[i+1] + u[i]))
    opti.subject_to( theta[i+1] == theta[i] + dt * (dtheta[i]))
    opti.subject_to( dtheta[i+1] == dtheta[i] + dt * (u[i] * cos(theta[i+1]) - sin(theta[i+1]) ))

opti.minimize( sum1(sin(theta)))

opti.solver("ipopt") #,p_opts, s_opts)
sol = opti.solve()
print(sol.value(x))
plt.plot(sol.value(x), label="x")
plt.plot(sol.value(u), label="u")
plt.plot(sol.value(theta), label="theta")
plt.legend()
plt.show()
'''
p = opti.parameter()
opti.set_value(p, 3)
'''

Very fast. Very impressive. Relatively readable code. I busted this out in like 15 minutes. IPopt solves the thing in the blink of an eye (about 0.05s self reported). Might be even faster if I warm start it with a good solution, as I would in online control (which may be feasible at this speed). Can add the initial condition as a parameter to the problem

I should try this on an openai gym.

Haskell has bindings to casadi.

http://hackage.haskell.org/package/dynobud

Thoughts on Faking Some of GADTs in Rust

I’m a guy who is somewhat familiar with Haskell who is trying to learn Rust. So I thought I’d try to replicate some cool Haskell functionality in Rust. I would love to hear comments, because I’m trying to learn. I have no sense of Rust aesthetics yet and in particular I have no idea how this interacts with the borrow system. What follows is a pretty rough brain dump

GADTs (Generalized algebraic data types) are an extension in Haskell that allows you to write constrained type signatures for your data constructors. They also change how the type checking of pattern matching is processed.

GADTs are sometimes described/faked by being built by making data types that hold equality/unification constraints. Equality constraints in Haskell like a ~ Int are fairly magical and the Rust compiler does not support them in an obvious way. Maybe this is the next project. Figure out how to fake ’em if one can. I don’t think this is promising though, because faking them will be a little wonky, and then GADTs are a little wonky on top of that. See https://docs.rs/refl/0.1.2/refl/ So we’ll go another (related) road.

This is roughly what GADTs look like in Haskell.

data MyGadt a where
  TBool :: Bool -> MyGadt Bool
  TInt :: Int -> MyGadt Int

And here is one style of encoding using smart constructors and a typeclass for elimination (pattern matching is replicated as a function that takes callbacks for the data held in the different cases). Regular functions can have a restricted type signature than the most general one their implementation implies. The reason to use a typeclass is so that we can write the eliminator as returning the same type that the GADT supplies. There isn’t an explicit equality constraint. A kind of Leibnitz equality

(forall f. f a -> f b)
is hiding in the eliminator. The Leibnitz equality can be used in place of (~) constraints at some manual cost. http://code.slipthrough.net/2016/08/10/approximating-gadts-in-purescript/

Click to access leibniz-equality.pdf


data MyGadt2 a = TBool2 Bool | TInt2 Int
-- smart constructors. Hide the original constructors
tInt :: Int -> MyGadt2 Int
tInt = TInt

tBool :: Bool -> MyGadt2 Bool
tBool = TBool

-- gadt eliminator
class MyGadtElim a where
   elimMyGadt :: forall f. MyGadt2 a -> (Bool -> f Bool) -> (Int -> f Int) -> f a

instance MyGadtElim Int where
   elimMyGadt (TInt2 x) fb fi = fi x
   elimMyGadt _ _ _ = error "How did TBool2 get type MyGadt2 Int?"
instance MyGadtElim Bool where
   elimMyGadt (TBool2 x) fb fi = fb x
 

The

forall f :: * -> *
is a problem for Rust. Rust does not have higher kinded types, although they can be faked to some degree. https://gist.github.com/CMCDragonkai/a5638f50c87d49f815b8 There are murmurs of Associated Type Constructors / GATs , whatever those are , that help ease the pain, but I’m pretty sure they are not implemented anywhere yet.

I’m going to do something related, a defunctionalization of the higher kinded types. We make an application trait, that will apply the given type function tag to the argument. What I’m doing is very similar to what happens in the singletons library, so we may be getting some things for free.

https://typesandkinds.wordpress.com/2013/04/01/defunctionalization-for-the-win/

trait App1 { type T;}

Then in order to define a new typelevel function rip out a quick tag type and an App impl.

struct F1 {}
impl App1 for F1{
    type T = Vec>;
}

It might be possible to sugar this up with a macro. It may also be possible to write typelevel functions in a point free style without defining new function tag names. The combinators Id, Comp, Par, Fst, Snd, Dup, Const are all reasonably definable and fairly clear for small functions. Also the combinator S if you want to talk SKI combinatory calculus, which is unfit for humans. https://en.wikipedia.org/wiki/SKI_combinator_calculus For currying, I used a number for how many arguments are left to be applied (I’m not sure I’ve been entirely consistent with these numbers). You need to do currying quite manually. It may be better to work with tuplized arguments

struct Vec1 {} // number is number of aspplications left.
impl App1 for Vec1{
    type T = Vec;
}

struct Id();

impl App1 for Id{
    type T = A; 
}

// partially applied const.
// struct Const()
// call this next one const1
struct Const1();
struct Const(PhantomData); // should be Const1
impl App1 for Const1 { // partial application
    type T = Const;
}

impl App1 for Const{
    type T = B; 
}



struct Fst {}
impl  App1<(A,B)> for Fst {
    type T = A;
} 

struct Snd {}
impl  App1<(A,B)> for Snd {
    type T = B;
} 


struct Dup{}
impl  App1 for Dup {
    type T = (A,A);
} 

struct Par2 {}
struct Par1 (PhantomData);
struct Par (PhantomData , PhantomData );
impl App1 for Par2 {
    type T = Par1;
}
impl App1 for Par1 {
    type T = Par;
}

impl App1<(X,Y)> for Par where F : App1, G : App1  {
    type T = (F::T, G::T);
}

// In order to curry, i think I'd have to define a name for every curried form.


// combinator calculus Const is K, Id is I, and here is S combinator. Yikes.
type I = Id;
type K = Const;

struct S3{}
struct S2(PhantomData);
struct S1(PhantomData, PhantomData);
// struct S(PhantomData, PhantomData, PhantomData);
impl  App1 for S1 where A : App1, B : App1, >::T : App1< >::T > {
    type T =  < >::T   as App1< >::T > >::T;
}  


struct Comp2();
struct Comp1 (PhantomData);
struct Comp (PhantomData, PhantomData);
impl App1 for Comp where G : App1, F : App1<>::T> {
    type T =  >::T>>::T;
}  

Anyway, the following is a translation of the above Haskell (well, I didn’t wrap an actual i64 or bool in there but I could have I think). You need to hide the actual constructors labeled INTERNAL in a user inaccessible module.

enum Gadt {
    TIntINTERNAL(PhantomData),
    TBoolINTERNAL(PhantomData)
    
}
// then build the specialzied constructors that 
fn TBool() -> Gadt{
    Gadt::TBoolINTERNAL(PhantomData)
}

fn TInt() -> Gadt{
    Gadt::TIntINTERNAL(PhantomData)
}

The smart constructors put the right type in the parameter spot

Then pattern matching is a custom trait per gadtified type. Is it possible to unify the different elimination traits that will come up into a single Elim trait? I’m 50-50 about whether this is possible. What we’re doing is a kind of fancy map_or_else if that helps you.

https://doc.rust-lang.org/std/option/enum.Option.html#method.map_or_else

trait GadtElim {
    type Inner;
fn gadtElim(&self, case1 : >::T   , case2 : >::T ) -> >::T  where 
         F : App1,
         F : App1,
         F : App1;
         

}


impl GadtElim for Gadt {
    type Inner = i64;
    fn gadtElim(&self, case1 : >::T   , case2 : >::T ) -> >::T  where 
         F : App1,
         F : App1,
         F : App1{
        match self{
            Gadt::TIntINTERNAL(PhantomData) => case2,
            Gadt::TBoolINTERNAL(PhantomData) => panic!("Somehow TBool has type Gadt")// Will never be reached though. god willing
        }
    }
}
impl GadtElim for Gadt {
    type Inner = bool;
    fn gadtElim(&self, case1 : >::T , case2 : >::T ) -> >::T  where 
         F : App1,
         F : App1,
         F : App1{
        match self{
            Gadt::TIntINTERNAL(PhantomData) => panic!("Somehow TInt has type Gadt"),
            Gadt::TBoolINTERNAL(PhantomData) => case1 // Will never be reached though. god willing
        }
    }
}

Usage. You have to explicitly pass the return type function to the eliminator. No inference is done for you. It’s like Coq’s match but worse. BTW the dbg! macro is the greatest thing on earth. Well done, Rust.

    let z = TInt().gadtElim::(true , 34);
    let z2 = TBool().gadtElim::(true , 34);

    dbg!(z);
    dbg!(z2);
    struct F7 {} // You need to do this. Kind of sucks. macroify? App!(F = Vec)
    impl App1 for F7 { type T = Vec; }
    dbg!(TInt().gadtElim::(vec!(true,false) , vec!(34,45,4,3,46)));

You can make helpers that don’t require explicit types to be given

fn gadtRec(x : impl GadtElim, case1 : A, case2 : A) -> A {
    x.gadtElim::>(case1 , case2)
}

One could also make an Eq a b type with Refl similarly. Then we need typelevel function tags that take two type parameter. Which, with currying or tupling, we may already have.

Questions:

Is this even good? Or is it a road of nightmares? Is this even emulating GADTs or am I just playing phantom type games?

We aren’t at full gadt. We don’t have existential types. Rust has some kind of existential story evolving (already there?), but the state of it is confusing to me. Something to play with. Higher rank functions would help?

Are overlapping typeclasses a problem in Rust?

Again, I have given nearly zero thought to borrowing and how it interacts with this. I’m a Rust n00b. I should think about it. Different eliminators based on whether you own or are borrowing?.

How much of singleton style dependent types do we get from this? It feels like we have already paid the cost of defunctionalizing. http://hackage.haskell.org/package/singletons

My current playground for this is at https://github.com/philzook58/typo

Edit: Rust is an ML with affine types and a syntax facelift roughly. So Ocaml is a good place to pump. Oleg Kiselyov had a fascinating approach that sort of smuggles through an Option type in an equality type using mutation. I wonder how well that would mesh with Rust. It seems obviously not thread safe unless you can make the mutation and matching atomic.

http://okmij.org/ftp/ML/GADT.txt

okmij.org/ftp/ML/first-class-modules/

https://blog.janestreet.com/more-expressive-gadt-encodings-via-first-class-modules/

Cvxpy and NetworkX Flow Problems

Networkx outputs scipy sparse incidence matrices

https://networkx.github.io/documentation/networkx-1.9/reference/generated/networkx.linalg.graphmatrix.incidence_matrix.html#networkx.linalg.graphmatrix.incidence_matrix

https://docs.scipy.org/doc/scipy/reference/sparse.html

Networkx also has it’s own flow solvers, but cvxpy gives you some interesting flexibility, like turning the problem mixed integer, quadratic terms, and other goodies. Plus it is very easy to get going as you’ll see.

So here’s a basic example of putting these two together. Very straightforward and cool.


import networkx as nx
import cvxpy as cvx
import matplotlib.pyplot as plt
import numpy as np
from scipy.sparse import lil_matrix

#graph is an networkx graph from somewhere

#print(edgedict)
nEdges = len(graph.edges)
nNodes = len(graph.nodes)

posflow = cvx.Variable(nEdges)
negflow = cvx.Variable(nEdges)

# split flow into positive and negative parts so we can talk about absolute value.
# Perhaps I should let cvxpy do it for me
constraints = [ 0 <= posflow,  0 <= negflow ]

absflow = posflow + negflow
flow = posflow - negflow

L = nx.incidence_matrix(graph, oriented=True )

source = np.zeros(nNodes) #lil_matrix(n_nodes)
# just some random source placement.
source[7] = 1
source[25] = -1

# cvxpy needs sparse matrices wrapped.
Lcvx = cvx.Constant(L)
#sourcecvx = cvx.Constant(source)

# flow conservation
constraints.append(Lcvx*flow == source)

# can put other funky inequality constraints on things.

objective = cvx.Minimize(cvx.sum(absflow)) 

print("building problem")
prob = cvx.Problem(objective, constraints)
print("starting solve")
prob.solve(solver=cvx.OSQP, verbose = True) #or try cvx.CBC, cvx.CVXOPT, cvx.GLPK, others
np.set_printoptions(threshold=np.inf)
print(absflow.value)


Here was a cool visual from a multi commodity flow problem (nx.draw_networkx_edges)

Nice, huh.

A Little Bloop on Typed Template Haskell

I’ve found looking at my statistics that short, to the point, no bull crap posts are the most read and probably most useful.

I’ve been tinkering around with typed template Haskell, which a cursory internet search doesn’t make obvious even exists. The first thing to come up is a GHC implementation wiki that seems like it might be stalled in some development branch. No. Typed template Haskell is already in GHC since GHC 7.8. And the first thing that comes up on a template haskell google search is the style of Template Haskell where you get deep into the guts of the syntax tree, which is much less fun.

Here’s a basic addition interpreter example.

{-# LANGUAGE TemplateHaskell -#}
import Language.Haskell.TH
import Language.Haskell.TH.Syntax

data Expr = Val Int | Add Expr Expr

eval' :: Expr -> Int
eval' (Val n) = n
eval' (Add e1 e2) = (eval' e1) + (eval' e2)

eval :: Expr -> TExpQ Int
eval (Val n) = [|| n ||]
eval (Add e1 e2) = [|| $$(eval e1) + $$(eval e2) ||]

The typed splice $$ takes a out of TExpQ a. The typed quote [|| ||] puts it in. I find that you tend to be able to follow the types to figure out what goes where. If you’re returning a TExpQ, you probably need to start a quote. If you need to recurse, you often need to splice. Etc.

You tend to have to put the actual use of the splice in a different file. GHC will complain otherwise.

ex1 :: Int
ex1 = $$(eval (Add (Val 1) (Val 1)))

At the top of your file put this to have the template haskell splices dumped into a file

{-# OPTIONS_GHC -ddump-splices -ddump-to-file #-}

Or have your package.yaml look something like this

executables:
  staged-exe:
    main:                Main.hs
    source-dirs:         app
    ghc-options:
    - -threaded
    - -rtsopts
    - -with-rtsopts=-N
    - -ddump-splices
    - -ddump-to-file
    dependencies:
    - staged

If you’re using stack, you need to dive into .stack/dist/x86/Cabal/build and then /src or into the executable folder something-exe-tmp/app to find .dump-splices files.

app/Main.hs:11:10-35: Splicing expression
    eval (Add (Val 1) (Val 1)) ======> (1 + 1)

Nice. I don’t know whether GHC might have optimized this thing down anyhow or not, but we have helped the process along.

Some more examples: unrolling the recursion on a power function (a classic)

-- ordinary version
power 0 x = 1
power n x = x * (power (n-1) x)

-- templated up
power' :: Int -> TExpQ (Int -> Int)
power' 0 = [|| const 1 ||]
power' n = [|| \x -> x * $$(power' (n-1)) x ||] 

You can unroll a fibonacci calculation

-- unrolled fib woth sharing
fib 0 = 1
fib 1 = 1
fib n = (fib (n-1)) + (fib (n-2))

-- we always need [||  ||] wheneve there is a Code
fib' :: Int -> TExpQ Int
fib' 0 = [|| 1 ||]
fib' 1 = [|| 1 ||] 
fib' n = [|| $$(fib' (n-1)) + $$(fib' (n-2)) ||]

This is blowing up in calculation though (no memoization, you silly head). We can implement sharing in the code using let expression (adapted from https://wiki.haskell.org/The_Fibonacci_sequence). Something to think about.

fib4 :: Int -> TExpQ Int
fib4 n = go n [|| ( 0, 1 ) ||]
                where
                  go :: Int -> TExpQ (Int, Int) -> TExpQ Int
                  go n z | n==0      = [|| let (a,b) = $$(z) in a ||]
                         | otherwise = go (n-1) [|| let (a,b) = $$(z) in (b, a + b) ||]

Tinkering around, you’ll occasionally find GHC can’t splice and quote certain things. This is related to cross stage persistence and lifting, which are confusing to me. You should look in the links below for more info. I hope I’ll eventually get a feel for it.

If you want to see more of my fiddling in context to figure out how to get stuff to compile here’s my github that I’m playing around in https://github.com/philzook58/staged-fun

Useful Link Dump:

Simon Peyton Jones has a very useful talk slides. Extremely useful https://www.cl.cam.ac.uk/events/metaprog2016/Template-Haskell-Aug16.pptx

Matthew Pickering’s post keyed me into that Typed Template Haskell is even there.

 http://mpickering.github.io/posts/2019-02-14-stage-3.html

MuniHac 2018: Keynote: Beautiful Template Haskell

https://markkarpov.com/tutorial/th.html Mark Karpov has a useful Template Haskell tutorial

Oleg Kiselyov: I’m still trying to unpack the many things that are going on here. Oleg’s site and papers are a very rich resource.

I don’t think that staged metaprogramming requires finally tagless style, as I first thought upon reading some of his articles. It is just very useful.

http://okmij.org/ftp/meta-programming/

http://okmij.org/ftp/meta-programming/tutorial/

Typed Template Haskell is still unsound (at least with the full power of the Q templating monad) http://okmij.org/ftp/meta-programming/#TExp-problem

You can write things in a finally tagless style such that you can write it once and have it interpreted as staged or in other ways. There is a way to also have a subclass of effects (Applicatives basically?) less powerful than Q that is sound.