Integration of values in a buffer with Delphi

Question

Following on from my question on differentiation:

I'm now looking at doing the integration. I can't quite get my head around this one. The situation is that I receive a buffer of data periodically that contains a number of values that are a fixed distance in time apart. I need to differentiate them. It is soo long since I did calculus at school ....

What I have come up with is this:

procedure IntegrateBuffer(ABuffer: TDoubleDynArray;
                                   var AOutBuffer: TDoubleDynArray;
                                   AVPS: integer);
const
  SumSum: double = 0.0;
  LastValue: double = NaN;
var
  i: integer;
  dt, aa, hl, hr: double;
begin
  // protect from divide by zero
  if (AVPS < 1) then exit;

  dt := 1 / AVPS;

  for i := 0 to high(ABuffer) do begin
    if (i = 0) then begin
      if (IsNaN(LastValue)) then begin
        hl := ABuffer[0];
        hr := ABuffer[0];
      end else begin
        hl := LastValue;
        hr := ABuffer[i];
      end;
    end else begin
      hl := ABuffer[i -1];
      hr := ABuffer[i];
    end;

    aa := 0.5 * dt * (hl + hr);
    SumSum := SumSum + aa;
    AOutBuffer[i] := SumSum;
  end;

  // remember the last value for next time
  LastValue := ABuffer[high(ABuffer)];
end;

I'm using the trapezium rule, hl and hr ar the left and right heights of the trapezium. dt is the base.

AVPS is values per second. A typical value for this would be between 10 and 100. The length of the buffers would typically be 500 to 1000 values.

I call the buffer time after time with new data which is continuous with the previous block of data, hence keeping the last value of the block for next time.

Is what I have done correct? ie, will it integrate the values properly?

Thank you.

What happened when you tested it? This is not a site for us to test your code for you. Clearly you are very competent. You must learn to write tests for yourself. Do you want me to write the exact same answer as I did last time? Also, not very keen on your use of writeable typed constants. Do not use them for this. — David Heffernan, Oct 04 '13 at 08:37
When I tested it, I didn't get values I was expecting. I went through the process of differentiating a series of numbers then put the results through this code and the numbers didn't match much. I don't know if I have the algorithm right before I code it. — user745323, Oct 04 '13 at 08:57
Your test is too complex. Test the function `y(x) = 1`. That should integrate to `b-a` where b and a are the integration limits. Then test `y(x) = x`. The indefinite integral of that is `1/2 x^2` so the definite integral is `(b^2 - a^2)/2`. Remember that trapezium is not exact so you'll need a small dt to get close the true values. Trapezium will be exact for these two functions though. Finally, if you've done testing, show it. Present an SSCCE with expected input and output. When testing always start with the simplest test you can think of. — David Heffernan, Oct 04 '13 at 09:01
You know, I have no clue how to do the tests you suggest. Furthermore, it is not the type of data I will be working with. And having done the tests and they come out wrong, as I suspect they will as I don't really understand the algorithm yet, which is why I'm asking for help, I wouldn't know what to do to correct the code. A series of 10 numbers similar to real life is a pretty simple test, and one that I understand. Thanks, though. — user745323, Oct 04 '13 at 09:28
Your algorithm looks ok. I might say something about the decision to use assignable typed constants and I might take issue with the fact that `SumSum` is going to continue to grow ad-infinitum as you continue to use this function (it will overflow eventually). That said, you say you tested this by differentiating a set of numbers, then re-integrating them. Numerical derivatives can be very tricky - they're hard to make precise and you almost always lose some information. I wouldn't expect even a good numerical derivative to re-integrate exactly all the time. David's tests are the answer. — J..., Oct 04 '13 at 10:12
@J... Actually `SumSum` won't overflow. Suppose that `aa` has value equal to `1`. Once `SumSum` reaches around `1e16`, `SumSum + aa == SumSum`. — David Heffernan, Oct 04 '13 at 10:20
@DavidHeffernan right... good point. Still, it will effectively break. — J..., Oct 04 '13 at 10:21
@J... Yes, that's true but that might not be a problem here. And it's a little tricky to fix. — David Heffernan, Oct 04 '13 at 10:23
@DavidHeffernan Perhaps we're talking about different problems... I was of the mind that even local writable typed constants are only initialized once (not on every method call). Is that not the case? — J..., Oct 04 '13 at 10:25
@J... Your understanding is correct. The idea is that this function is called repeatedly with new data and extends the integration starting where it left off at the end of the previous call. So writeable typed constant works. But it's lame. Doesn't allow reset. It's a very poor solution. — David Heffernan, Oct 04 '13 at 10:34
@DavidHeffernan Yes, that was my only point - that this function would eventually run out of usable lifespan as it continued to be used. — J..., Oct 04 '13 at 10:36

David Heffernan · Answer 1 · 2013-10-04T10:37:38.387

Looks like you need some help with testing the code. Here, as discussed in comments, is a very simple test.

{$APPTYPE CONSOLE}

uses
  SysUtils, Math;

type
  TDoubleDynArray = array of Double;

var
  SumSum: double;
  LastValue: double;

procedure Clear;
begin
  SumSum := 0.0;
  LastValue := NaN;
end;

procedure IntegrateBuffer(
  ABuffer: TDoubleDynArray;
  var AOutBuffer: TDoubleDynArray;
  AVPS: integer
);
var
  i: integer;
  dt, aa, hl, hr: double;
begin
  // protect from divide by zero
  if (AVPS < 1) then exit;

  dt := 1 / AVPS;

  for i := 0 to high(ABuffer) do begin
    if (i = 0) then begin
      if (IsNaN(LastValue)) then begin
        hl := ABuffer[0];
        hr := ABuffer[0];
      end else begin
        hl := LastValue;
        hr := ABuffer[i];
      end;
    end else begin
      hl := ABuffer[i -1];
      hr := ABuffer[i];
    end;

    aa := 0.5 * dt * (hl + hr);
    SumSum := SumSum + aa;
    AOutBuffer[i] := SumSum;
  end;

  // remember the last value for next time
  LastValue := ABuffer[high(ABuffer)];
end;

var
  Buffer: TDoubleDynArray;
  OutBuffer: TDoubleDynArray;

begin
  // test y = 1 for a single call, expected output = 1, actual output = 2
  Clear;
  Buffer := TDoubleDynArray.Create(1.0, 1.0);
  SetLength(OutBuffer, Length(Buffer));
  IntegrateBuffer(Buffer, OutBuffer, 1);
  Writeln(OutBuffer[high(OutBuffer)]);

  Readln;
end.

I'm integrating the function y(x) = 1 over the range [0..1]. So, the expected output is 1. But the actual output is 2.

So, what's wrong? You can work it out in the debugger, but it's easy enough to see by inspecting the code. You are summing a triangle on the very first sample. When IsNaN(LastValue) is true then you should not make a contribution to the integral. At that point you've not covered any distance on the x axis.

So to fix the code, let's try this:

....
if (IsNaN(LastValue)) then begin
  hl := 0.0;//no contribution to sum
  hr := 0.0;
end else begin
  hl := LastValue;
  hr := ABuffer[i];
end;
....

That fixes the problem.

Now let's extend the test a little and test y(x) = x:

// test y = x, expected output = 12.5
Clear;
Buffer := TDoubleDynArray.Create(0.0, 1.0, 2.0, 3.0, 4.0, 5.0);
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Writeln(OutBuffer[high(OutBuffer)]);

So, that looks good.

OK, what about multiple calls:

// test y = x for multiple calls, expected output = 18
Clear;
Buffer := TDoubleDynArray.Create(0.0, 1.0);
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Buffer := TDoubleDynArray.Create(2.0, 3.0, 4.0, 5.0, 6.0);
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Writeln(OutBuffer[high(OutBuffer)]);

And how about one value at a time?

// test y = x for multiple calls, one value at a time, expected 0.5
Clear;
Buffer := TDoubleDynArray.Create(0.0);
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Buffer := TDoubleDynArray.Create(1.0);
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Writeln(OutBuffer[high(OutBuffer)]);

What about passing an empty array?

// test y = x for multiple calls, some empty arrays, expected 0.5
Clear;
Buffer := TDoubleDynArray.Create(0.0);
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Buffer := nil;
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Buffer := TDoubleDynArray.Create(1.0);
SetLength(OutBuffer, Length(Buffer));
IntegrateBuffer(Buffer, OutBuffer, 1);
Writeln(OutBuffer[high(OutBuffer)]);

Uh, oh, access violation. Better protect that by simply skipping the function at the start if the buffer is empty:

if (AVPS < 1) then exit;
if (Length(ABuffer) = 0) then exit;

OK, now that last test passes

Hopefully you get the idea now. I've just used noddy Writeln based testing but that does not scale. Get yourself a unit test framework (I recommend DUnitX) and build proper test cases. This will also force you to factor your code so that it is well designed. One of the often unexpected benefits of making code testable is that it usually results in the design of the interface being improved.

For your next question, I request that you supply an SSCCE with the test code! ;-)

Some comments on the code:

Pass dynamic arrays by const or by var. In your case you want to pass the input buffer by const.
Don't use writeable typed constants. Use either parameters, or some other more sane state management.

Again, as I said in the previous question, write tests to prove code, as well as checking it by eye. The key to writing tests is to start with the very simplest thing you can possibly think of. Something so simple that you know for 100% sure the answer. Then, once you get that to work, expand the testing to more complex cases.

Thank you David. That is very helpful. My reasoning behind setting hl and hr to ABuffer[0] first time through is that I figured that as the real life data won't be starting at zero, a good place to start would be by using the first value in the array. I still think this is a reasonable idea, but it does fail your tests. I'll muse some more on that. As for writeable constants, what your issue with them? I use them as though they are persistent local variables. They are not needed outside of the method, needn't be stored outside the method. They encourage encapsulation. Thanks again. — user745323, Oct 04 '13 at 10:58
@user745323 They make your code untestable for a start. That's because there's no way to re-initialise those variables. If you want to maintain state from call to call make the function a method of an object, or pass in the state as a parameter. — David Heffernan, Oct 04 '13 at 11:01
As for what to do with the first value, it's up to you. Your code assumes that the integration starts from time `-dt`, and that the signal value was constant between time `-dt` and time `0`. My change makes assumes that integration starts from `0`. — David Heffernan, Oct 04 '13 at 11:02
Thanks. I think my assumption is closer to real life, but it will probably have negligible effect in practice. — user745323, Oct 04 '13 at 11:06
David, just to confirm, you tested only the last value return, are the other values in the array valid, too? Thanks. (I wish I was better at maths) — user745323, Oct 04 '13 at 11:10
As for your assumption on the first value it comes down to this. If the first sample you get is the start point, my code is right. If the first sample you get is second value, and for some reason you miss the first value, your code is right. — David Heffernan, Oct 04 '13 at 12:06
@user745323 - Wishes are for fairy tales. If you want to be better at maths I suggest you find a cheap used text and start refreshing. You say it's been a while since you studied calculus but at least you've studied it and already have something of a foot in the door. It may not have seemed relevant at the time but now that it is evident that it is important for what you do I'd say it is well worth the time and effort to put in some practice and sharpen those skills again. It really can't be stated strongly enough how important solid maths are to progressing as a good programmer. — J..., Oct 04 '13 at 12:40
@David, regarding the assumption - I'm converting from acceleration to velocity. The thing I'm measuring may already have some acceleration by the time I start getting the values. I can't assume the thing is starting from rest. Its not really a case of missing a value. — user745323, Oct 04 '13 at 12:57
You are assuming it is starting from rest. You assume it is at rest at time `t=-dt`. — David Heffernan, Oct 04 '13 at 13:03

Integration of values in a buffer with Delphi

1 Answers1