Skip to content

Commit

Permalink
Merge branch 'master' of github.com:psg-mit/micromarshall
Browse files Browse the repository at this point in the history
  • Loading branch information
bmsherman committed Jul 7, 2023
2 parents 51c167b + ad43130 commit f9d00cf
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 79 deletions.
50 changes: 25 additions & 25 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# smooth

This is the artifact for the paper "λS: Computable semantics for differentiable programming with higher-order functions and datatypes".
This is the artifact for the paper "<img src="https://render.githubusercontent.com/render/math?math=\lambda_S">: Computable semantics for differentiable programming with higher-order functions and datatypes".

This repository contains the implementation of <img src="https://render.githubusercontent.com/render/math?math=\lambda_S"> as an embedded language within Haskell.
We name this library _smooth_.
Expand Down Expand Up @@ -52,8 +52,8 @@ let gradient (f : ℜ^2 → ℜ) (x : ℜ^2) : ℜ^2 =
```
! camera assumed to be at the origin
let raytrace (s : Surface (ℜ^2)) (lightPos : ℜ^2) (rayDirection : ℜ^2) : ℜ =
let t_star = firstRoot (λ t : ℜ ⇒ s (scale t rayDirection)) in
let y = scale t_star rayDirection in let normal : ℜ^2 = - gradient s y in
let tStar = firstRoot (λ t : ℜ ⇒ s (scale t rayDirection)) in
let y = scale tStar rayDirection in let normal : ℜ^2 = - gradient s y in
let lightToSurf = y - lightPos in
max 0 (dot (normalize normal) (normalize lightToSurf))
/ (norm2 y * norm2 lightToSurf)
Expand Down Expand Up @@ -90,12 +90,12 @@ raytrace :: VectorSpace g => ((DReal :* DReal) :=> DReal) g ->
(DReal :* DReal) g ->
(DReal :* DReal) g -> DReal g
raytrace s lightPos u =
let t = firstRoot (ArrD (\wk t -> dmap wk s # (scale t (dmap wk u)))) in
let y = scale t u in
let tStar = firstRoot (ArrD (\wk t -> dmap wk s # (scale t (dmap wk u)))) in
let y = scale tStar u in
let normal = gradient s y in
let lightVector = lightPos `sub` y in
max 0 (dot (normalize normal) (normalize lightVector))
/ (norm2 y * norm2 lightVector)
let lightToSurf = lightPos `sub` y in
max 0 (dot (normalize normal) (normalize lightToSurf))
/ (norm2 y * norm2 lightToSurf)
where
(x0 :* x1) `sub` (y0 :* y1) = (x0 - y0) :* (x1 - y1)

Expand All @@ -122,14 +122,14 @@ runRaytraceExampleDeriv = atPrec 1e-3 raytraceExampleDeriv
#### Section 2.1
Code in <img src="https://render.githubusercontent.com/render/math?math=\lambda_S">:
```
let t_star y = firstRoot (λ t : ℜ ⇒ 1 - y^2 - (t - 1)^2)
let tStar y = firstRoot (λ t : ℜ ⇒ 1 - y^2 - (t - 1)^2)
```
```
deriv t_star y = - deriv (λ y0 : ℜ ⇒ f (t_star y) y0) y /
deriv (λ t : ℜ ⇒ f t y) (t_star y)
deriv tStar y = - deriv (λ y0 : ℜ ⇒ f (tStar y) y0) y /
deriv (λ t : ℜ ⇒ f t y) (tStar y)
```
```
deriv t_star y = - y / (t_star y - 1):
deriv tStar y = - y / (tStar y - 1):
```
Implementation in _smooth_:
```haskell
Expand Down Expand Up @@ -286,27 +286,27 @@ secondDerivVarianceLinearChange =
Code in <img src="https://render.githubusercontent.com/render/math?math=\lambda_S">:

```
eps=1e-3> hausdorffDist d_R2 l_shape (quarterCircle 0)
eps=1e-3> hausdorffDist distR2 lShape (quarterCircle 0)
[0.4138, 0.4145]
```
```
eps=1e-1> deriv (λ y : ℜ ⇒ hausdorffDist d_R2 l_shape (quarterCircle y)) 0
eps=1e-1> deriv (λ y : ℜ ⇒ hausdorffDist distR2 lShape (quarterCircle y)) 0
[-0.752, -0.664]
```

Implementation in _smooth_:
```haskell
-- Section 7.3: Hausdorff distance between quarter-circle and L-shape.
quarter_circle :: VectorSpace g => DReal g -> Maximizer (DReal :* DReal) g
quarter_circle y0 = M.map (ArrD (\wk theta -> let y0' = dmap wk y0 in
(cos (pi / 2 * theta)) :* (y0' + sin (pi / 2 * theta)))) M.unit_interval
quarterCircle :: VectorSpace g => DReal g -> Maximizer (DReal :* DReal) g
quarterCircle y0 = M.map (ArrD (\wk theta -> let y0' = dmap wk y0 in
(cos (pi / 2 * theta)) :* (y0' + sin (pi / 2 * theta)))) M.unitInterval

quarter_square_perim :: VectorSpace g => Maximizer (DReal :* DReal) g
quarter_square_perim = M.union (M.map (ArrD (\_ x -> x :* 1)) M.unit_interval)
(M.map (ArrD (\_ y -> 1 :* y)) M.unit_interval)
lShape :: VectorSpace g => Maximizer (DReal :* DReal) g
lShape = M.union (M.map (ArrD (\_ x -> x :* 1)) M.unitInterval)
(M.map (ArrD (\_ y -> 1 :* y)) M.unitInterval)

hausdorffDistCircleL :: DReal ()
hausdorffDistCircleL = hausdorffDist d_R2 quarter_square_perim (quarter_circle 0)
hausdorffDistCircleL = hausdorffDist distR2 lShape (quarterCircle 0)

-- Time: 7 seconds
-- Result: [0.41384921849465670653300003702, 0.41440539111235744884709353399]
Expand All @@ -317,10 +317,10 @@ runHausdorffDistCircleL = atPrec 0.001 hausdorffDistCircleL
-- Section 7.3: derivative of Hausdorff distance between quarter-circle and L-shape wrt. translation.
hausdorffDistTranslatedQuarterCircle :: DReal ()
hausdorffDistTranslatedQuarterCircle =
deriv (ArrD (\_ y -> hausdorffDist d_R2 quarter_square_perim (quarter_circle y))) 0
deriv (ArrD (\_ y -> hausdorffDist distR2 lShape (quarterCircle y))) 0

-- Time: 52 seconds
-- Result: [-0.752, -0.664]
-- Result: [-0.7515973045396820224886373089321421844, -0.6641561255883687886832219076605117364]
runHausdorffDistTranslatedQuarterCircle :: Real
runHausdorffDistTranslatedQuarterCircle = atPrec 0.1 hausdorffDistTranslatedQuarterCircle
```
Expand All @@ -336,7 +336,7 @@ We provide both a virtual machine image that can be directly downloaded and run

### Virtual machine image

We provide a virtual machine image with Ubuntu 20.04 with all of the dependencies preloaded, packaged together with this README file as SmoothVM.ova. You can import the .ova file into hypervisor software (e.g., virtualbox).
We provide a [virtual machine image](https://drive.google.com/file/d/1xx30qckdGbW1Sukn8G_f1_tJZEKKoorq/view?usp=sharing) with Ubuntu 20.04 with all of the dependencies preloaded, packaged together with this README file as SmoothVM.ova. You can import the .ova file into hypervisor software (e.g., virtualbox).
Note that the README file in the source code within that VM is out of date; please prefer using this README for evaluation.

Once the virtual machine is loaded, you can sign in to the user `lambda-s` with the password `lambda-s`. Open a terminal and use the command `cd smooth` to access the repository. View the examples from the paper with `vim src/SmoothLang.hs`.
Expand Down Expand Up @@ -414,4 +414,4 @@ Objects of the category **HAD** in the paper directly correspond to Haskell func
The real type is the type family `DReal` defined in `src/Real.hs`.
Products are defined by the type family `:*` in `src/Experimental/PSh.hs`.
Exponentials are defined by the type family `:=>` in `src/FwdPSh.hs`.
Various derived types with **HAD** are defined in files in the directory `src/Types/`.
Various derived types with **HAD** are defined in files in the directory `src/Types/`.
30 changes: 15 additions & 15 deletions src/SmoothLang.hs
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ import Types.Real
import Types.Integral (Integral, mean, variance, uniform)
import FwdMode (getValue, VectorSpace)
import Rounded (ofString, RoundDir( Down ))
import Types.Maximizer (hausdorffDist, d_R2, Maximizer)
import Types.Maximizer (hausdorffDist, distR2, Maximizer)
import qualified Types.Maximizer as M

dRealtoReal :: DReal () -> CPoint Real
Expand Down Expand Up @@ -51,12 +51,12 @@ raytrace :: VectorSpace g => ((DReal :* DReal) :=> DReal) g ->
(DReal :* DReal) g ->
(DReal :* DReal) g -> DReal g
raytrace s lightPos u =
let t = firstRoot (ArrD (\wk t -> dmap wk s # (scale t (dmap wk u)))) in
let y = scale t u in
let tStar = firstRoot (ArrD (\wk t -> dmap wk s # (scale t (dmap wk u)))) in
let y = scale tStar u in
let normal = gradient s y in
let lightVector = lightPos `sub` y in
max 0 (dot (normalize normal) (normalize lightVector))
/ (norm2 y * norm2 lightVector)
let lightToSurf = lightPos `sub` y in
max 0 (dot (normalize normal) (normalize lightToSurf))
/ (norm2 y * norm2 lightToSurf)
where
(x0 :* x1) `sub` (y0 :* y1) = (x0 - y0) :* (x1 - y1)

Expand Down Expand Up @@ -163,16 +163,16 @@ secondDerivVarianceLinearChange =


-- Section 7.3: Hausdorff distance between quarter-circle and L-shape.
quarter_circle :: VectorSpace g => DReal g -> Maximizer (DReal :* DReal) g
quarter_circle y0 = M.map (ArrD (\wk theta -> let y0' = dmap wk y0 in
(cos (pi / 2 * theta)) :* (y0' + sin (pi / 2 * theta)))) M.unit_interval
quarterCircle :: VectorSpace g => DReal g -> Maximizer (DReal :* DReal) g
quarterCircle y0 = M.map (ArrD (\wk theta -> let y0' = dmap wk y0 in
(cos (pi / 2 * theta)) :* (y0' + sin (pi / 2 * theta)))) M.unitInterval

quarter_square_perim :: VectorSpace g => Maximizer (DReal :* DReal) g
quarter_square_perim = M.union (M.map (ArrD (\_ x -> x :* 1)) M.unit_interval)
(M.map (ArrD (\_ y -> 1 :* y)) M.unit_interval)
lShape :: VectorSpace g => Maximizer (DReal :* DReal) g
lShape = M.union (M.map (ArrD (\_ x -> x :* 1)) M.unitInterval)
(M.map (ArrD (\_ y -> 1 :* y)) M.unitInterval)

hausdorffDistCircleL :: DReal ()
hausdorffDistCircleL = hausdorffDist d_R2 quarter_square_perim (quarter_circle 0)
hausdorffDistCircleL = hausdorffDist distR2 lShape (quarterCircle 0)

-- Time: 7 seconds
-- Result: [0.41384921849465670653300003702, 0.41440539111235744884709353399]
Expand All @@ -183,9 +183,9 @@ runHausdorffDistCircleL = atPrec 0.001 hausdorffDistCircleL
-- Section 7.3: derivative of Hausdorff distance between quarter-circle and L-shape wrt. translation.
hausdorffDistTranslatedQuarterCircle :: DReal ()
hausdorffDistTranslatedQuarterCircle =
deriv (ArrD (\_ y -> hausdorffDist d_R2 quarter_square_perim (quarter_circle y))) 0
deriv (ArrD (\_ y -> hausdorffDist distR2 lShape (quarterCircle y))) 0

-- Time: 52 seconds
-- Result: [-0.752, -0.664]
-- Result: [-0.7515973045396820224886373089321421844, -0.6641561255883687886832219076605117364]
runHausdorffDistTranslatedQuarterCircle :: Real
runHausdorffDistTranslatedQuarterCircle = atPrec 0.1 hausdorffDistTranslatedQuarterCircle
12 changes: 6 additions & 6 deletions src/Types/KShape.hs
Original file line number Diff line number Diff line change
Expand Up @@ -86,19 +86,19 @@ quarter_disk_variable :: VectorSpace g => DReal g -> KShape (DReal :* DReal) g
quarter_disk_variable r = intersect unit_square $
ArrD $ \wk (x :* y) -> x^2 + y^2 < (dmap wk r)^2

d_R2 :: ((DReal :* DReal) :* (DReal :* DReal) :=> DReal) g
d_R2 = ArrD $ \_ ((x :* y) :* (x' :* y')) -> (x - x')^2 + (y - y')^2
distR2 :: ((DReal :* DReal) :* (DReal :* DReal) :=> DReal) g
distR2 = ArrD $ \_ ((x :* y) :* (x' :* y')) -> (x - x')^2 + (y - y')^2

d_R1 :: (DReal :* DReal :=> DReal) g
d_R1 = ArrD $ \_ (x :* x') -> (x - x')^2
distR1 :: (DReal :* DReal :=> DReal) g
distR1 = ArrD $ \_ (x :* x') -> (x - x')^2

exampleHausdorffDist :: DReal ()
exampleHausdorffDist = hausdorffDist d_R2 unit_square quarter_disk
exampleHausdorffDist = hausdorffDist distR2 unit_square quarter_disk

-- Seems to be converging, but slowly. Derivatives
-- must not be working, because it's only progressing through bisection.
exampleHausdorffDist2 :: DReal ()
exampleHausdorffDist2 = hausdorffDist d_R1 unit_interval unit_interval
exampleHausdorffDist2 = hausdorffDist distR1 unit_interval unit_interval

xPlusY :: ((DReal :* DReal) :=> DReal) g
xPlusY = ArrD (\_ (x :* y) -> x + y)
Expand Down
66 changes: 33 additions & 33 deletions src/Types/Maximizer.hs
Original file line number Diff line number Diff line change
Expand Up @@ -51,34 +51,34 @@ separationDist :: VectorSpace g => PShD a =>
separationDist d k k' =
inf k' (ArrD (\wk x' -> inf (dmap wk k ) (ArrD (\wk' x -> dmap (wk @. wk') d # (x :* dmap wk' x')))))

unit_interval :: VectorSpace g => Maximizer DReal g
unit_interval = ArrD $ \_ -> max01
unitInterval :: VectorSpace g => Maximizer DReal g
unitInterval = ArrD $ \_ -> max01

unit_square :: VectorSpace g => Maximizer (DReal :* DReal) g
unit_square = product unit_interval unit_interval
unit_square = product unitInterval unitInterval

quarter_disk :: VectorSpace g => Maximizer (DReal :* DReal) g
quarter_disk = quarter_disk_variable 1

quarter_disk_variable :: VectorSpace g => DReal g -> Maximizer (DReal :* DReal) g
quarter_disk_variable r = compactUnion unit_interval $ ArrD (\wk z ->
compactUnion unit_interval $ ArrD (\wk' theta ->
quarter_disk_variable r = compactUnion unitInterval $ ArrD (\wk z ->
compactUnion unitInterval $ ArrD (\wk' theta ->
let r' = dmap (wk @. wk') r in
let z' = dmap wk' z in
point ((r' * z' * cos ((pi / 2) * theta)) :* (r' * z' * sin ((pi / 2) * theta)))))

d_R2 :: ((DReal :* DReal) :* (DReal :* DReal) :=> DReal) g
d_R2 = ArrD $ \_ ((x :* y) :* (x' :* y')) -> sqrt ((x - x')^2 + (y - y')^2)
distR2 :: ((DReal :* DReal) :* (DReal :* DReal) :=> DReal) g
distR2 = ArrD $ \_ ((x :* y) :* (x' :* y')) -> sqrt ((x - x')^2 + (y - y')^2)

d_R1 :: (DReal :* DReal :=> DReal) g
d_R1 = ArrD $ \_ (x :* x') -> (x - x')^2
distR1 :: (DReal :* DReal :=> DReal) g
distR1 = ArrD $ \_ (x :* x') -> (x - x')^2

circle :: VectorSpace g => DReal g -> Maximizer (DReal :* DReal) g
circle r = map (ArrD (\wk theta -> let r' = dmap wk r in
(r' * cos (2 * pi * theta)) :* (r' * sin (2 * pi * theta)))) unit_interval
(r' * cos (2 * pi * theta)) :* (r' * sin (2 * pi * theta)))) unitInterval

unit_square_perim :: VectorSpace g => Maximizer (DReal :* DReal) g
unit_square_perim = foldr1 union [ map f unit_interval | f <- fs ]
unit_square_perim = foldr1 union [ map f unitInterval | f <- fs ]
where
fs :: VectorSpace g => [(DReal :=> (DReal :* DReal)) g]
fs = [ ArrD (\_ x -> (2 * x - 1) :* (-1))
Expand All @@ -87,14 +87,14 @@ unit_square_perim = foldr1 union [ map f unit_interval | f <- fs ]
, ArrD (\_ y -> 1 :* (2 * y - 1))]

exampleHausdorffDist :: DReal ()
exampleHausdorffDist = hausdorffDist d_R2 unit_square quarter_disk
exampleHausdorffDist = hausdorffDist distR2 unit_square quarter_disk

exampleHausdorffDist2 :: DReal ()
exampleHausdorffDist2 = hausdorffDist d_R1 unit_interval unit_interval
exampleHausdorffDist2 = hausdorffDist distR1 unitInterval unitInterval

-- Does this run?
exampleHausdorffDistDeriv :: DReal ()
exampleHausdorffDistDeriv = deriv (ArrD (\_ r -> hausdorffDist d_R2 unit_square (quarter_disk_variable r))) 1
exampleHausdorffDistDeriv = deriv (ArrD (\_ r -> hausdorffDist distR2 unit_square (quarter_disk_variable r))) 1

xPlusY :: ((DReal :* DReal) :=> DReal) g
xPlusY = ArrD (\_ (x :* y) -> x + y)
Expand All @@ -106,25 +106,25 @@ exampleMaximizationDeriv :: DReal ()
exampleMaximizationDeriv = deriv (ArrD (\_ r -> sup (quarter_disk_variable r) xPlusY)) 1

simplerMaximization :: DReal ()
simplerMaximization = supremum (map (ArrD (\_ x -> 0.5 * x)) unit_interval)
simplerMaximization = supremum (map (ArrD (\_ x -> 0.5 * x)) unitInterval)

simpleDerivTest :: DReal ()
simpleDerivTest = deriv (ArrD (\_ c -> supremum ((map (ArrD (\wk x -> dmap wk c * x))) unit_interval))) 1.0
simpleDerivTest = deriv (ArrD (\_ c -> supremum ((map (ArrD (\wk x -> dmap wk c * x))) unitInterval))) 1.0


car :: VectorSpace g => DReal g -> Maximizer (DReal :* DReal) g
car y = map (ArrD (\wk theta -> (cos (pi * theta) :* (sin (pi * theta) + dmap wk y)))) unit_interval
car y = map (ArrD (\wk theta -> (cos (pi * theta) :* (sin (pi * theta) + dmap wk y)))) unitInterval

obstacle :: VectorSpace g => Maximizer (DReal :* DReal) g
obstacle = union
(map (ArrD (\_ x -> ((2 + x) :* 2))) unit_interval)
(map (ArrD (\_ y -> (2 :* (2 + y)))) unit_interval)
(map (ArrD (\_ x -> ((2 + x) :* 2))) unitInterval)
(map (ArrD (\_ y -> (2 :* (2 + y)))) unitInterval)

carClearance :: VectorSpace g => DReal g
carClearance = separationDist d_R2 (car 0) obstacle
carClearance = separationDist distR2 (car 0) obstacle

derivCarClearance :: VectorSpace g => DReal g
derivCarClearance = deriv (ArrD (\_ y -> separationDist d_R2 (car y) obstacle)) 0
derivCarClearance = deriv (ArrD (\_ y -> separationDist distR2 (car y) obstacle)) 0

cvx2 :: DReal g -> (DReal :* DReal) g -> (DReal :* DReal) g -> (DReal :* DReal) g
cvx2 c (x0 :* x1) (y0 :* y1) =
Expand All @@ -133,8 +133,8 @@ cvx2 c (x0 :* x1) (y0 :* y1) =

triangleR2 :: VectorSpace g => (DReal :* DReal) g -> (DReal :* DReal) g -> (DReal :* DReal) g -> Maximizer (DReal :* DReal) g
triangleR2 x y z =
compactUnion unit_interval (ArrD (\wk a ->
compactUnion unit_interval (ArrD (\wk' b ->
compactUnion unitInterval (ArrD (\wk a ->
compactUnion unitInterval (ArrD (\wk' b ->
let wk'' = wk @. wk' in
point (cvx2 (dmap wk' a) (dmap wk'' x) (cvx2 b (dmap wk'' y) (dmap wk'' z)))))
))
Expand All @@ -152,26 +152,26 @@ obstacleCvx = convexHullR2 (point (2 :* 2) `union` point (2 :* 4) `union` point

carCvx :: VectorSpace g => Maximizer (DReal :* DReal) g
carCvx =
compactUnion unit_interval (ArrD (\_ r ->
map (ArrD (\wk theta -> (dmap wk r * cos (pi * theta)) :* (dmap wk r * sin (pi * theta)))) unit_interval
compactUnion unitInterval (ArrD (\_ r ->
map (ArrD (\wk theta -> (dmap wk r * cos (pi * theta)) :* (dmap wk r * sin (pi * theta)))) unitInterval
))
`union`
compactUnion unit_interval (ArrD (\_ x ->
compactUnion unit_interval (ArrD (\wk y ->
compactUnion unitInterval (ArrD (\_ x ->
compactUnion unitInterval (ArrD (\wk y ->
point ((-dmap wk x) :* (-3 * y))
))
))

carClearanceCvx :: VectorSpace g => DReal g
carClearanceCvx = separationDist d_R2 (car 0) obstacleCvx
carClearanceCvx = separationDist distR2 (car 0) obstacleCvx

simpleObstacle :: VectorSpace g => Maximizer (DReal :* DReal) g
simpleObstacle = union
(map (ArrD (\_ x -> ((1 + x) :* 1))) unit_interval)
(map (ArrD (\_ y -> (1 :* (1 + y)))) unit_interval)
(map (ArrD (\_ x -> ((1 + x) :* 1))) unitInterval)
(map (ArrD (\_ y -> (1 :* (1 + y)))) unitInterval)

d_R2_squared :: ((DReal :* DReal) :* (DReal :* DReal) :=> DReal) g
d_R2_squared = ArrD $ \_ ((x :* y) :* (x' :* y')) -> (x - x')^2 + (y - y')^2
distR2_squared :: ((DReal :* DReal) :* (DReal :* DReal) :=> DReal) g
distR2_squared = ArrD $ \_ ((x :* y) :* (x' :* y')) -> (x - x')^2 + (y - y')^2

simplef :: (DReal :=> DReal) g
simplef = ArrD $ \_ y0 -> separationDist d_R2_squared (point (0 :* y0)) simpleObstacle
simplef = ArrD $ \_ y0 -> separationDist distR2_squared (point (0 :* y0)) simpleObstacle

0 comments on commit f9d00cf

Please sign in to comment.