Zulip Chat Archive

Stream: lean4

Topic: Making a RungeKutta 4

Miguel Raz Guzmán Macedo (Jun 17 2021 at 05:44):

Hello all!
I'm trying to dive again into Lean4.
I was inspired by this post to make a functional RungeKutta4 function https://www.johndcook.com/blog/2016/06/02/ode-solver-as-a-functional-fold/

So far, I have this:

import Data
def rk (t, y) t' : (t', y):= do
  let h  := t' - t
  let k1 := f t y
  let k2 := f (t + 0.5*h) (y + 0.5*h*k1)
  let k3 := f (t + 0.5*h) (y + 0.5*h*k2)
  let k4 := f (t + 1.0*h) (y + 1.0*h*k3)
  return (t', y + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0)

def f (t y : Float) : Float := (t^2 - y*2) * sin y

#eval f 0.1 0.2

#eval foldl rk (0.0, -1.0) #[0.01, 0.03, 0.04, 0.06]

My problems are the following:

I don't know how to import sin
I don't know how to use foldl
I don't know how to use the unit testing facilities

Miguel Raz Guzmán Macedo (Jun 17 2021 at 05:44):

Please ping if you are answering, thanks!

Mario Carneiro (Jun 17 2021 at 06:20):

Does that import Data line actually work for you? I'm getting unknown package Data but that might be my setup

Mario Carneiro (Jun 17 2021 at 06:24):

This works:

def rk (f : Float → Float → Float) : Float × Float → Float → Float × Float
| (t, y), t' =>
  let h  := t' - t
  let k1 := f t y
  let k2 := f (t + 0.5*h) (y + 0.5*h*k1)
  let k3 := f (t + 0.5*h) (y + 0.5*h*k2)
  let k4 := f (t + 1.0*h) (y + 1.0*h*k3)
  return (t', y + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0)

def f (t y : Float) : Float := (t^2 - y*2) * y.sin

#eval f 0.1 0.2

#eval Array.foldl (rk f) (init := (0.0, -1.0)) #[0.01, 0.03, 0.04, 0.06]

Mario Carneiro (Jun 17 2021 at 06:25):

if you want to open sin instead of using projection notation, you can do this

open Float (sin)
open Array (foldl)

def f (t y : Float) : Float := (t^2 - y*2) * sin y

#eval f 0.1 0.2

#eval foldl (rk f) (init := (0.0, -1.0)) #[0.01, 0.03, 0.04, 0.06]

Mario Carneiro (Jun 17 2021 at 06:26):

foldl takes its arguments in the opposite order, so you can either use named arguments like init := here or else swap the arguments

Mario Carneiro (Jun 17 2021 at 06:31):

and in case you are wondering why you are getting -1.1099 instead of -0.9527 like the blog post, it's because your definition of f uses y*2 instead of y^2

Miguel Raz Guzmán Macedo (Jun 17 2021 at 15:04):

Thanks a ton @Mario Carneiro , super helpful as always.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 15:05):

Next question: do I _need_ to specify the type signature of rk that explicitly?

Mario Carneiro (Jun 17 2021 at 20:13):

You still have to say that f is a function because lean does not try to guess if a function is dependent or not, but other than that you can do without all the types

def rk (f : _ → _ → _) := fun (t, y) t' =>
  let h  := t' - t
  let k1 := f t y
  let k2 := f (t + 0.5*h) (y + 0.5*h*k1)
  let k3 := f (t + 0.5*h) (y + 0.5*h*k2)
  let k4 := f (t + 1.0*h) (y + 1.0*h*k3)
  (t', y + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0)

Mario Carneiro (Jun 17 2021 at 20:15):

the return is not necessary and is bad for type inference

Mario Carneiro (Jun 17 2021 at 20:15):

@Miguel Raz Guzmán Macedo

Miguel Raz Guzmán Macedo (Jun 17 2021 at 20:20):

Oooh, those are very useful tips, thanks a lot! I didn't know about return.

Mario Carneiro (Jun 17 2021 at 20:25):

Normally you would use return in a do block, although it's not necessary at the last statement

Mario Carneiro (Jun 17 2021 at 20:26):

If you are just creating lets then you don't need do

Miguel Raz Guzmán Macedo (Jun 17 2021 at 20:26):

Ok, Gotcha.

Mario Carneiro (Jun 17 2021 at 20:27):

lean 4 has where syntax as well, but it doesn't seem to work with local variables introduced by pattern matching

Mario Carneiro (Jun 17 2021 at 20:30):

def rk (f : _ → _ → _) t y t' :=
  (t', y + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0)
where
  h  := t' - t
  k1 := f t y
  k2 := f (t + 0.5*h) (y + 0.5*h*k1)
  k3 := f (t + 0.5*h) (y + 0.5*h*k2)
  k4 := f (t + 1.0*h) (y + 1.0*h*k3)

Miguel Raz Guzmán Macedo (Jun 17 2021 at 20:30):

@Mario Carneiro nice, I had only seen the where for structs and types, so I didn't know how to use it.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 20:51):

Hmmm - I can't get the final line to work - here is all my code:

open Float (sin)
open Array (foldl)
def rk (f : _ → _ → _) t y t' :=
  (t', y + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0)
where
  h  := t' - t
  k1 := f t y
  k2 := f (t + 0.5*h) (y + 0.5*h*k1)
  k3 := f (t + 0.5*h) (y + 0.5*h*k2)
  k4 := f (t + 1.0*h) (y + 1.0*h*k3)

def f (t y : Float) : Float := (t^2 - y^2) * sin y

#eval f 0.1 0.2

#eval foldl rk (init := (0.0, -1.0)) #[0.01, 0.03, 0.04, 0.06]

Mario Carneiro (Jun 17 2021 at 21:02):

the last version changes the type of rk, so you would have to use uncurry (rk f) instead of rk in the eval line

Miguel Raz Guzmán Macedo (Jun 17 2021 at 21:24):

#eval foldl uncurry (rk f) (init := (0.0, -1.0)) #[0.01, 0.03, 0.04, 0.06]

Am I missing an export?

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:33):

The manual really needs some updating.

Mario Carneiro (Jun 17 2021 at 23:33):

(uncurry (rk f)) needs to be in parentheses since otherwise you are passing uncurry and (rk f) separately to foldl; also uncurry probably doesn't exist in lean 4 yet

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:33):

I keep wanting a Lean4 by example guide

Mario Carneiro (Jun 17 2021 at 23:33):

that's called the lean 4 source code ATM :upside_down:

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:34):

Even looking at the source code is not too helpful - how do you navigate when you want to see what a function does?

Mario Carneiro (Jun 17 2021 at 23:34):

ctrl-click

Mario Carneiro (Jun 17 2021 at 23:34):

Use vscode with the lean 4 extension to look at the lean 4 sources

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:35):

I'm used to searching for unit tests, but I don't think Lean4 has a good system for that.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:35):

@Mario Carneiro I have Lean3 installed in VSCode and Lean4 on Emacs, and that's after a few days of battle with the nix stuff.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:35):

I don't really get how to keep the installations separate.

Mario Carneiro (Jun 17 2021 at 23:35):

There are unit tests in lean 4 but not enough to be usable as an example guide

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:36):

thanks for your patience, sorry for being a newbie :(

Mario Carneiro (Jun 17 2021 at 23:36):

I have the lean 4 extension disabled by default and enabled on specific workspaces (like the lean4 repo)

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:36):

Oh that's nice.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:37):

I just want a little package that I can make little examples for.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:37):

Like exercism.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:37):

And start to benchmark different solutions.

Mario Carneiro (Jun 17 2021 at 23:37):

I ran into some issues with lean files getting misinterpreted as lean4 if I don't selectively enable the extension

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:37):

Understandable.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:38):

@Mario Carneiro A bit besides the point, but how likely is Lean5 to completely change Lean4 semantics? Do you find much transition pain from Lean3 to Lean4?

Mario Carneiro (Jun 17 2021 at 23:38):

I have no idea whether lean 5 will exist or what it will do

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:39):

@Mario Carneiro hoo boy.

Mario Carneiro (Jun 17 2021 at 23:39):

there is a lot of transition pain from lean 3 to lean 4, but working on fixing that is literally in my job description at the moment

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:39):

@Mario Carneiro Oh neat - where at?

Mario Carneiro (Jun 17 2021 at 23:40):

microsoft research

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:40):

Oh neat, so I guess there's at least 3 people working on Lean4 now? Kudos!

Mario Carneiro (Jun 17 2021 at 23:41):

heh, I'm not on the core team, but my spirit is

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:41):

Neat.

Mario Carneiro (Jun 17 2021 at 23:42):

I'm working with Daniel Selsam on the mathport tool

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:42):

@Mario Carneiro Neat!

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:42):

Wait, which extension is the good Lean4 extension?

Mario Carneiro (Jun 17 2021 at 23:43):

mhuisi.lean4 I think

Mario Carneiro (Jun 17 2021 at 23:44):

no wait, it's leanprover.lean4 now

Mario Carneiro (Jun 17 2021 at 23:46):

By the way, it's easy to write uncurry yourself:

def uncurry (f : α → β → γ) : α × β → γ | (a, b) => f a b
#eval Array.foldl (uncurry (rk f)) (init := (0.0, -1.0)) #[0.01, 0.03, 0.04, 0.06]

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:47):

Oh neat.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:47):

Ugh, I'm stuck in elan hell again.

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:51):

Alright, not worth my sanity to fight this package manager today. thanks though!

Miguel Raz Guzmán Macedo (Jun 17 2021 at 23:54):

That uncurry doesn't work btw

Mario Carneiro (Jun 17 2021 at 23:56):

full example:

def rk (f : _ → _ → _) t y t' :=
  (t', y + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0)
where
  h  := t' - t
  k1 := f t y
  k2 := f (t + 0.5*h) (y + 0.5*h*k1)
  k3 := f (t + 0.5*h) (y + 0.5*h*k2)
  k4 := f (t + 1.0*h) (y + 1.0*h*k3)

def f (t y : Float) : Float := (t^2 - y^2) * y.sin

#eval f 0.1 0.2

def uncurry (f : α → β → γ) : α × β → γ | (a, b) => f a b

#eval Array.foldl (uncurry (rk f)) (init := (0.0, -1.0)) #[0.01, 0.03, 0.04, 0.06]