bouncingball - "drag" ?

[miguel-negrao]

Letting the example "bouncingball" run for a while it behaves as if it had drag, the height that the ball reaches gets lower and lower until it doesn't move. And then it starts going down until disappears from the window. Both the Yampa and Bearriver versions behave the same. Is this due to loss of precision in the integration ? Is this a known issue ?

[turion]

How long does it take? I'm not completely surprised because it could be a rounding error, but I didn't know before.

[ivanperez-keera]

It is a know issue. You can check the numbers being produced, and yes, I think it's a rounding error.

[miguel-negrao]

It takes more or less 1 minute until the issue arises. It feels a bit strange that one of the main examples of the library cannot run a simulation for more than a couple of seconds. I tried changing to Double but that didn't help. How do people creating games using the library deal with this issue ?

[walseb]

I would add a conditional to stop it once its reached 0.

[miguel-negrao]

Can you clarify what do you mean ? I believe the simulation is supposed to be a fully elastic collision, without losses of any kind, so it should go on jumping to the same initial height forever, right ?

[walseb]

I see I thought the loss of force was intentional but not the negative force. Sorry for the confusion.

[ivanperez-keera]

I think one way to go about it is to detect the actual point of collision. Right now, the point / time where the collision with 0 does happen is probably not being calculated properly, so the bounce is not fully elastic. There's a small error, which is smaller the faster you sample, but it's there.

[miguel-negrao]

I see. The paper "Functional Reactive Programming, Refactored" touches on this and provides a solution with CCDGameMonad. Then I would suggest changing this example either to the most simple version of collision detection that works correctly in Yampa/Bearriver, or change to something without collision such as a pendulum.

[ivanperez-keera]

This example is useful because it's valid in both Dunai and Yampa. The example degenerates to the ideal meaning if the sampling rate is infinitely high. (I can also tell you that a game I built in the past used this, so it can be used --with small workarounds.)

Yampa cannot accept the custom monad, so that solution does not solve all problems.

But I still agree that it may be confusing, and it (rightfully) prompted this question. So perhaps 1) a comment in the file is in order, and 2) a small example with code that does not exhibit this bug could also be added separately.

[miguel-negrao]

I gave it a go at estimating the velocity at the real point of impact (in between ticks) and then estimating the final velocity for that tick after the collision, but I wasn't successful.

sf = bouncingBall ((300.0 :: Double), 0.0)

bouncingBall (p0,v0) =
  switch (proc (_) -> do
            let a = -99.8
            v <- (v0 +) ^<< integral -< a
            p <- (p0 +) ^<< integral -< v
            bounce <- edge               -< (p <= 0 && v < 0)
            pPrev <- iPre p0 -< p
            vPrev <- iPre v0 -< v
            t <- time -< ()
            tPrev <- iPre 0 -< t
            let
              qEquation (a, b, c) = (x1, x2)
                where
                  x1 = e + sqrt d / (2 * a)
                  x2 = e - sqrt d / (2 * a)
                  d = b * b - 4 * a * c
                  e = - b / (2 * a) 
              dt = (t - tPrev)
                                                                -- estimating under velocity vPrev the time delta at which collision occurs
              dtAtCollision = snd $ qEquation (a, vPrev, pPrev) -- p == 0 <=> pPrev + dt v == 0 <=> pPrev + dt (vPrev + dt a) == 0 <=> a dt² + vPrev dt + pPrev == 0
                                                                -- solve for dt 
              vAtCollision =  vPrev + (dtAtCollision * a)
              dtAfterCollision = dt - dtAtCollision
              vAfterCollision = (- vAtCollision) + (dtAfterCollision * a)
              pAfterCollision = trace (printf "vAtCollision: %.6f" vAtCollision) $ dtAfterCollision * vAfterCollision
            returnA -< ((p,v), bounce `tag` (pAfterCollision, vAfterCollision))
         )
         bouncingBall

vAtCollision: -244.205843
vAtCollision: -243.208437
vAtCollision: -242.211111
vAtCollision: -241.213854
vAtCollision: -240.216490
vAtCollision: -239.219360

[ivanperez-keera]

I suspect you can demonstrate the error by:

• Fixing the time step (e.g., 1s).

• Setting up the experiment such that the ball does not collide exactly at a time step (that is, the ball is already in the ground when the collision takes place).

In theory, the vertical velocity right before bouncing and right after bouncing should be the negation of each other. Also, the position should be the exact same.

If I am right, the experiment will show that they are not.

[ScottFreeCode]

I just saw this at random, but, I believe I found a while back that the integration formula is wrong. The formula used is: change in distance equals time passed times the speed at the current time (which is itself calculated from the total time passed).

freeFall :: Monad m => BallVel -- ^ The start velocity
  -> BehaviourF m UTCTime () Ball
freeFall v0 = arr (const (0, 0, -9.81))
  >>> integralFrom v0 >>> integral

This is a minor physics error. It should be: change in distance equals time passed times the average speed over that time duration (i.e., since acceleration is constant and thus difference in speed over time is linear, halfway between the speed at the current time and the speed at the previous time).

freeFall :: Monad m => BallVel -- ^ The start velocity
  -> BehaviourF m UTCTime () Ball
freeFall v0 = arr (const (0, 0, -9.81))
  >>> integralFrom v0 >>> averageLast v0 >>> integral

averageLast :: Monad m => BallVel -> BehaviourF m UTCTime BallVel BallVel
averageLast init = arr (^/2) >>> arr id &&& iPre (zeroVector ^+^ init ^/ 2) >>> arr (uncurry (^+^))

For clarity (and because I'm not certain I implemented the vectorspace version right), a scalar version of averageLast:

averageLast init = arr (/2) >>> arr id &&& iPre (init / 2) >>> arr (uncurry (+))

I believe this bug went back to the original Dunai paper, though the code there is just slightly different in implementation. My code is from the Rhine paper example, which I was adapting to test a particular extension. (I need to circle back around on that in general but it hasn't been high priority for a couple reasons I'll get back to.)

[miguel-negrao]

I've tested @ScottFreeCode 's solution and it fixes the issue. See my PR. I also fixed the time, but it is not related with this, if you just pass the delta as 0.01, and use averageLast, that still fixes the issue. I've added a black "floor" and visually the collision seems to be pretty accurate. I guess to make this really correct we would have to study a bit numerical integration of differential equations, there are many methods.

[ivanperez-keera]

I think this goes back to the original Yampa work. It's great that you identified this issue and a fix.

I'll work on including these as I wait for answers on other issues.

[ivanperez-keera]

Ok, just to make my life simpler, I've re-written the expression as:

averageLast v0 = (arr id &&& iPre v0) >>^ (\(x, y) -> (x + y) / 2)

Making a single change to dunai-examples/bouncingball/BouncingBall.hs so that it reads:

fallingBall p0 v0 = proc () -> do                                               
  v <- (v0 +) ^<< integral -< (-99.8)                                           
  p <- (p0 +) ^<< integral <<< averageLast v0 -< v                              
  returnA -< (p, v)                                                             

averageLast v0 = (arr id &&& iPre v0) >>^ (\(x, y) -> (x + y) / 2) 

seems to exhibit the right behavior.

However, when I read this expression, what I understand is that it is taking the speed to integrate in between sampling points, as opposed to at the end of the sampling period. I have a bit of a problem with this because, if my understanding is correct, then this error stems from the fact that the sampling frequency is limited. If sampling frequency was infinitely high, the error would not manifest. Or, as the time delta goes to zero (in the limit), the behavior converges to the ideal.

You have to be aware that there is sampling to determine how to write it, and the FRP abstraction breaks at this point.

So I'm just wondering if the original equation was wrong in the ideal continuous-time semantics, or if it was wrong because we are sampling (in practice).

[miguel-negrao]

I believe that by using averageLast we are actually numerically integrating the velocity using the trapezoidal rule instead of the rectangle rule. I think in general the trapezoidal rule will be more accurate. Actually, is there a reason why integral is not implemented using the trapezoidal rule ? I don't see any drawbacks in doing that, and it would then be more accurate, I think.

In general any correct numerical integration method will converge to the right result when the interval size converges to zero, but some methods converge faster, and might be therefore preferable, although computational complexity has to be taken into account also.

[ivanperez-keera]

It has to do with meaning more than implementation.

The selling point of FRP is that people should be able to program without thinking of sampling. This is a case where that abstraction breaks down.

Ideally, we'd put the solution in the definition of integral itself, but we can't, because the trapezoidal rule will only be numerically accurate (up to floating point precision) if the value integrated changes linearly. Otherwise, it may break (it would be a coincidence if it didn't).

One can fix a specific example, but I'm interested in thinking about how to best tell this story, and if anything can be done so that users of bearriver/Yampa can program without being aware of sampling.

Maybe it makes sense to define an integrateLinear that is the same as integral in the limit, but numerically more accurate for linearly changing inputs in practice. We'd be characterizing the input as a whole based on the shape of the derivative, which is something that exists for all time (it's a continuous-time notion). Or perhaps there's a better solution.

Something to think about.

[miguel-negrao]

Ideally, we'd put the solution in the definition of integral itself, but we can't, because the trapezoidal rule will only be numerically accurate (up to floating point precision) if the value integrated changes linearly. Otherwise, it may break (it would be a coincidence if it didn't).

Hum, I don't quite understand this. The trapezoidal rule is valid for any integrable function, same as the rectangle rule, not just for linear functions. I believe that in general for functions in C^2 (with second derivative which is continuous) the trapezoidal rule will converge faster with O(N^-2) while the rectangle rule converges with only O(N^-1) where N is the number of intervals. For functions which are not C^2 the analysis would have to be done differently.

For doing physics without worrying about sampling I think you would need to use symbolic calculation, i.e., actually solving the equations mathematically. Off course, this can only be done in limited cases. A programming language that does that is for instance "Mathematica". Otherwise I think we will always need to worry somewhat about sampling and the methods used when approximating integrals or derivatives, as we are using numerical methods underneath.

[ivanperez-keera]

Sure, the rule is valid in the sense that you can use it, but like the rectangle rule, the result is just an approximation of the integral.

I still think it could make sense to transport the thinking to the input signal first. Ideally, yes, the programming language would have knowledge of the particulars of the signal being integrated and do so symbolically. In absence, the user has to do it.

Perhaps giving the user integration functions that are characterized by properties of input they integrate could be one way to go. Or maybe we really have to provide several orders or Runge-Kuttas. Idk.

[miguel-negrao]

In the comment in your commit you have "Alternative to 'integrate' that is numerically more accurate when integrating linearly changing inputs.". I don't see why it matters if the function is linear or not, but maybe I'm missing something. If it is linear then the trapezoid rule should give the exact result and not just an approximation. If it is not linear but is C^2 it still converges faster with smaller dt to the right result than the rectangle rule.

[ivanperez-keera]

Whether these approximation methods will be more accurate or not will depend on the method, the function being integrated, and the sampling rate.

In this case I'm not speaking about covergence. I assume the user cannot lower the sampling rate at will.

I don't think it's true that trapezoid will be more accurate than rectangle rule regardless of the sampling rate (it will be more accurate for many sampling rates, especially if they are uniform, but not necessarily all). For some functions you could probably pick sampling rates where the result of integral is closer to the ideal than the result of integralLinear.

Additionally, wrt. to the phrase "numerically more accurate" and not just saying "exact": The output is still subject to floating-point rounding errors. This is implemented as a sum, so it's numerically more accurate, but not perfectly accurate. Not that this is different from the lack of precision at a specific point. It's not the case in general for floating point numbers that x + y using FP arithmetic is the same as the floating-point representation of the "ideal" number x+y. The former will carry higher imprecision (the one in the representation of x, the one in the representation of y, and the imprecision of FP addition).

[miguel-negrao]

I got curious about this and built some code to test the accuracy of rectangle vs trapezoid method. The code is available here. The error is calculated using the mean square error. I tested with a couple of functions with easy to calculate antiderivatives. The error values are low enough to give some confidence that I calculated correctly the analytical integral. The integral is calculated using the fundamental theorem of calculus integral f(x) dx from 0 to t = F(t) - F(0) where F is one primitive or antiderivative of f. The sampling happens with fixed delta. All the intervals of the functions where chosen to get enough change going in that time period.

Results:

numPoints: 20
                   Description:  Trap bttr dur   MSE Rect   MSE Trap      ratio
                      f(x) = 1:       True  10    0.000e0    0.000e0        NaN
                      f(x) = x:       True  10    2.135e0    0.000e0   Infinity
                   f(x) = x**2:       True  10    1.400e2   5.932e-2    2.359e3
                   f(x) = x**3:       True  10    1.065e4    8.402e0    1.267e3
                   f(x) = x**6:       True  10    6.890e9    1.752e7    3.932e2
    f(x) = -x**3+10 x**2+4 x+1:       True  15    1.869e4    4.929e1    3.792e2
                 f(x) = sin(x):       True  13   4.856e-2   1.567e-3    3.100e1
                    f(x) = e^x:       True   2   2.232e-2   5.997e-6    3.722e3
          (e**(-x/10))* cos(x):       True   8   4.358e-2   4.144e-5    1.052e3

numPoints: 100
                   Description:  Trap bttr dur   MSE Rect   MSE Trap      ratio
                      f(x) = 1:       True  10    0.000e0    0.000e0        NaN
                      f(x) = x:       True  10   8.375e-2  3.272e-27   2.560e25
                   f(x) = x**2:       True  10    5.117e0   9.306e-5    5.499e4
                   f(x) = x**3:       True  10    3.704e2   1.269e-2    2.919e4
                   f(x) = x**6:       True  10    2.074e8    2.376e4    8.731e3
    f(x) = -x**3+10 x**2+4 x+1:       True  15    5.790e2   7.165e-2    8.081e3
                 f(x) = sin(x):       True  13   1.957e-3   2.573e-6    7.605e2
                    f(x) = e^x:       True   2   8.189e-4   9.039e-9    9.060e4
          (e**(-x/10))* cos(x):       True   8   1.774e-3   6.799e-8    2.610e4

numPoints: 1000
                   Description:  Trap bttr dur   MSE Rect   MSE Trap      ratio
                      f(x) = 1:       True  10    0.000e0    0.000e0        NaN
                      f(x) = x:       True  10   8.338e-4  1.946e-25   4.283e21
                   f(x) = x**2:       True  10   5.012e-2   9.264e-9    5.410e6
                   f(x) = x**3:       True  10    3.585e0   1.252e-6    2.863e6
                   f(x) = x**6:       True  10    1.938e6    2.283e0    8.488e5
    f(x) = -x**3+10 x**2+4 x+1:       True  15    5.448e0   7.005e-6    7.777e5
                 f(x) = sin(x):       True  13   1.972e-5  2.595e-10    7.599e4
                    f(x) = e^x:       True   2   8.028e-6  8.914e-13    9.006e6
          (e**(-x/10))* cos(x):       True   8   1.783e-5  6.845e-12    2.605e6

numPoints: 10000
                   Description:  Trap bttr dur   MSE Rect   MSE Trap      ratio
                      f(x) = 1:       True  10    0.000e0    0.000e0        NaN
                      f(x) = x:       True  10   8.334e-6  7.977e-24   1.045e18
                   f(x) = x**2:       True  10   5.001e-4  9.260e-13    5.401e8
                   f(x) = x**3:       True  10   3.573e-2  1.250e-10    2.858e8
                   f(x) = x**6:       True  10    1.925e4   2.274e-4    8.464e7
    f(x) = -x**3+10 x**2+4 x+1:       True  15   5.414e-2  6.989e-10    7.746e7
                 f(x) = sin(x):       True  13   1.974e-7  2.597e-14    7.599e6
                    f(x) = e^x:       True   2   8.012e-8  8.902e-17    9.000e8
          (e**(-x/10))* cos(x):       True   8   1.784e-7  6.849e-16    2.604e8

ratio is MSE rect / MSE trap . The error with the rectangle method is between 3x larger and up to 25 orders of magnitude larger ! The error with rectangle is larger even when dropping the number of points to 20. which is quite to low for 2 periods of the sine function.

The functions used for testing here are some of the most common functions, sines, polynomials, exponentials. Some more complicated functions could be used, for instance getting the solution to a differential equations, calculating the derivative of that and integrating numerically. I can try to get some examples of that sort.

But I have to say that just conceptually I have a hard time thinking of a function where the rectangle rule would perform better. Lets assume the function is continuous. If the sampling is so coarse that the function is wiggling up and down between two time steps then, any of the two methods will be almost like a random number, any slightly different sampling will give a completely different result. If the sampling is fine enough that the between two time steps the the function is monotonic (either goes up or goes down) with some type of curve then clearly the trapezoid will approach the integral better then a rectangle. The only exception is the constant function where both methods give a perfect result.

In the end I believe what matters is which of the methods gives the best result in most cases, in terms of typical functions used in dunai or FRP.