yampa: improve accuracy of integral

[idontgetoutmuch]

Luckily I found this https://github.com/ivanperez-keera/dunai/issues?q=is%3Aissue+is%3Aclosed. Without it the bouncing ball does not bounce for long.

EDIT (by Ivan Perez): This issue has been transferred from the dunai repo to the Yampa repo.

[ivanperez-keera]

Which specific issue were you trying to point to?

Is this a problem only in dunai or also in Yampa?

[ivanperez-keera]

@idontgetoutmuch Please let us know the details. I'll close this issue in a couple of days otherwise. Thanks!

[idontgetoutmuch]

With

fallingBall :: Monad m => Pos -> Vel -> SF m () (Pos, Vel)
fallingBall p0 v0 = proc () -> do
  let g = -9.81
  v <- arr (\x -> v0 + x) <<< integral -< g
  p <- arr (\y -> p0 + y) <<< integral -< v
  returnA -< (p,v)

bouncingBall :: Monad m => Double -> Double -> SF m () (Double, Double)
bouncingBall p0 v0 =
  switch (fallingBall p0 v0 >>> (arr id &&& hitFloor))
         (\(p,v) -> bouncingBall p (-v))

hitFloor :: Monad m => SF m (Double,Double) (Event (Double,Double))
hitFloor = arr $ \(p,v) ->
  if p < 0 && v < 0 then Event (p,v) else noEvent

I get

But with

fallingBall :: Monad m => Pos -> Vel -> SF m () (Pos, Vel)
fallingBall p0 v0 = proc () -> do
  let g = -9.81
  v <- arr (\x -> v0 + x) <<< integral -< g
  p <- arr (\y -> p0 + y) <<< integralLinear v0 -< v
  returnA -< (p,v)

integralLinear :: Monad m => Double -> SF m Double Double
integralLinear initial = average >>> integral
  where
    average = (arr id &&& iPre initial) >>^ (\(x, y) -> (x ^+^ y) ^/ 2)

I get what I expect

Full code is here: https://github.com/idontgetoutmuch/phd-jonathan-thaler/commits/literate-haskell-for-review. In particular, it should be possible to run this: https://github.com/idontgetoutmuch/phd-jonathan-thaler/blob/literate-haskell-for-review/SIRABM.lhs#L262.

[ivanperez-keera]

Thanks!

Based on past simulations, I think this is a problem with Yampa. If we can identify the problem there, we should fix it in Yampa and then just make bearriver adopt the same solution as Yampa (or one that is extensionally equivalent).

[idontgetoutmuch]

This (lifted from https://github.com/miguel-negrao/test-integral/blob/main/app/Main.hs) works but I don't know how you could apply this to Yampa (thanks @miguel-negrao).

integral :: (Monad m, VectorSpace a s) => SF m a a
integral = integralFrom zeroVector

integralFrom :: (Monad m, VectorSpace a s) => a -> SF m a a
integralFrom a0 =
  zeroVector --> proc a -> do
    dt <- constM ask -< ()
    aPrev <- iPre zeroVector -< a
    accumulateWith (^+^) a0 -< (realToFrac dt / 2) *^ (a ^+^ aPrev)

[idontgetoutmuch]

Maybe this?

integral :: forall a s . (Fractional s, VectorSpace a s) => SF a a
integral = SF {sfTF = tf0}
  where
    tf0 a0 = (integralAux igrl0 a0, igrl0)

    igrl0  = zeroVector

    integralAux :: a -> a -> SF' a a
    integralAux igrl a_prev = SF' tf -- True
      where
        tf :: DTime -> a -> (SF' a a, a)
        tf dt a = (integralAux igrl' a, igrl')
          where
            igrl' = igrl ^+^ (realToFrac dt / 2) *^ (a ^+^ a_prev)

[miguel-negrao]

Isn't this a duplicate of ivanperez-keera/dunai#268 ?

[idontgetoutmuch]

Isn't this a duplicate of ivanperez-keera/dunai#268 ?

Probably

[idontgetoutmuch]

I made the change and with a test program it works.

Here's my test

fallingBall :: Double -> Double -> SF () (Double, Double)
fallingBall p0 v0 = proc () -> do
  let g = -9.81
  v <- arr (\x -> v0 + x) <<< integral -< g
  p <- arr (\y -> p0 + y) <<< integral -< v
  returnA -< (p,v)

bouncingBall :: Double -> Double -> SF () (Double, Double)
bouncingBall p0 v0 =
  switch (fallingBall p0 v0 >>> (arr id &&& hitFloor))
         (\(p,v) -> bouncingBall p (-v))

hitFloor :: SF (Double, Double) (Event (Double, Double))
hitFloor = arr $ \(p,v) ->
  if p < 0 && v < 0 then Event (p,v) else noEvent

Here's the problem

and after the fix

But if I run the yampa tests then they just hang :-(

I have put a PR here so you can try it yourself: #262

[idontgetoutmuch]

I don't know what to do about tests that just hang. I don't understand yampa sufficiently to debug it.

[ivanperez-keera]

This issue pertains to Yampa. I'm transferring it there.

[ivanperez-keera]

There are some small errors (the height going over 15 or under 0). What happens if you run this for a really long time? Do those always cancel out or do they accumulate?

[chansey97]

@idontgetoutmuch Sorry, this may be a bit off-topic. Could you please tell me how you generate these signal diagrams with Yampa. Very thanks.

[idontgetoutmuch]

@chansey97 Full code is here: https://github.com/idontgetoutmuch/phd-jonathan-thaler/commits/literate-haskell-for-review. In particular, it should be possible to run this: https://github.com/idontgetoutmuch/phd-jonathan-thaler/blob/literate-haskell-for-review/SIRABM.lhs#L262.

[ivanperez-keera]

I'd like to move this issue forward.

Would someone be willing to help by putting together a test that measures the difference between the old integral and the new integral for several input functions with known integrals, and with varying sampling rates?

I'm thinking of something that will measure the distance between each implementation and the known integral, for example, measure arr cos >>> integral vs arr sin, calculate the average and max errors, and do it for both versions of the integral, for several functions (not just cos), and with varying sampling rates, and show that way that the new implementation is just more robust.

[vitalibarozzi]

I wrote a quick naive version of it. Was this more or less what you meant?

{-# LANGUAGE RankNTypes #-}
module Lib
    ( someFunc
    ) where

import FRP.Yampa
import FRP.Yampa.Integration
import qualified FRP.Yampa as Yampa
import Control.Monad

someFunc :: IO ()
someFunc = do

    let input  = replicate 100 (0.01, Nothing)

    let sinSF  = Yampa.time >>> arr sin
    let sinSFy = embed sinSF ((), input)

    putStrLn "\nOld version"

    let cosSF  = Yampa.time >>> arr cos >>> integral
    let cosSFy = embed cosSF  ((), input)
    let diff1  = map (\(a,b) -> abs $ a - b) (zip cosSFy sinSFy)

    putStrLn $ "Average error: " <> show (errAvg diff1)
    putStrLn $ "Max error    : " <> show (errMax diff1)

    putStrLn "\nNew version"

    let cosSF2  = Yampa.time >>> arr cos >>> integral_v2
    let cosSFy2 = embed cosSF2 ((), input)
    let diff2   = map (\(a,b) -> abs $ a - b) (zip cosSFy2 sinSFy)

    putStrLn $ "Average error: " <> show (errAvg diff2)
    putStrLn $ "Max error    : " <> show (errMax diff2)

  where

    errAvg xs = sum xs / realToFrac (length xs)
    errMax xs = foldr max 0 xs

Old version
Average error: 7.923829007273358e-4
Max error : 2.2914762007648637e-3

New version
Average error: 3.827574186695313e-6
Max error : 7.0122698941910144e-6

As far as I can tell the new version is better for most if not all cases.

[idontgetoutmuch]

I don't have time to think about it further at the moment but why not make the integration method a function that can be passed in (with a default)? I am thinking here about symplectic methods for differential equations. Also we ought to be able to estimate the error for given methods anyway (although a test is also nice). There will be a performance / accuracy pay-off of course.