sampling endless streaming timeserie into an array

[davinci]

Hi Fortran Folks,

I'm experimenting with the most efficient way to follow and sample high-frequency streaming timeseries into an array of fixed size.
This could be more of an "algorithms and datastructure" issue (Donald Knuth like), but as I expect, a typical problem for Fortran timeseries processing and hardware processor functionality too;

My Fortran app should process 10.000 arrays of size 10.000 each (or one type containing all arrays).

Every split of a second a new data-row is sampled (e.g. from an instrument or from stock exchanges) and inserted into the next free row. After the last row is reached, the oldest row is deleted and the new data can be inserted somehow. And so on, first out, last in. In time, the ORDERED array content looks kindlike 123, 234, 345...

---[1] I could shift the whole array before each insert at the end;

A(1:n-1) = A(2:n)
A(n) = newdata

---[2] Or implement some pointered linked list construction without array shift, replacing the oldest rows;

Time A1,2,3 P1,2,3

1 1 . . . . .
2 1 2 . . 1 .
3 1 2 3 . 1 2
4 4 2 3 3 . 2
5 4 5 3 3 1 .
6 4 5 6 . 1 2
7 7 5 6 3 . 2
8 7 8 6 3 1 .
9 7 8 9 . 1 2

In this way I have only one insertion in A and two in P each time.

---[3] Or use some smarter (faster) allocate/deallocate/pointer construction

---[4] Or some implicite function unknown to me, which is compiled to an extremely fast low level processing, using special (Intel) processor tasks and/or buffered I/O in the background ? Buffered I/O could be a requirement anyway.

---[5] Parallel HPF on quadcores (of course always an extra option for later).

Extreme efficiency is essential, and therefore as few instructions as possible.
I suppose there is a lot of experience in the Fortran world about algorithms and benchmarking of this kind of time series sampling, but couldn't find much yet.
Ideas and experiences are very welcome.
I'll post the results.

Thanks, Clemens de Leeuw (DaVinci)

[Jugoslav_Dujic]

An additional question is, of course, what are you going to do later with so obtained data structure? I'm asking because you may invent a blazing fast algorithm for sampling but ultimately end up having data in a format not easy (or slow) to convert.

In any case, option [1] is out of the question. Out of others (I'm not an expert on parallel execution), two schemes come to my mind. One is "plain vanilla" (nSamples,nSeries) matrix with reusable columns, with an additional "first row" indicator which changes as the new series arrive:

integer, parameter:: nSamples = 10000, nSeries=10000

real:: A(nSamples, nSeries)

integer:: ifirstcol

ifirstcol = 1

do while(something interrupts)

 call ReadData(A(:,ifirstcol))

 ifirstcol = mod(ifirstcol+1, nSeries) + 1

end do

In this way, your most recent data, in order, are A(:,ifirstcol), A(:,ifirstcol+1),..., A(:,nSeries), A(:,1), ...,A(:ifirstcol-1). IOW, you simply rewrite the oldest column each time, located at ifirstcol.

I'd venture a guess that this may end up as being the most efficient way. An alternative approach could be either with an array(nSeries) or a linked list of derived TYPEs, each containing an array of (nSamples) elements, or a circularly linked list of those same derived TYPEs. In linked list approach, you simply change the head pointer to the next one each time you receive the data (much the same as ifirstcol changes):

integer, parameter:: nSamples = 10000, nSeries=10000

type t_mydata

 type(t_mydata), pointer:: next

 real:: data(nSamples)

end type t_mydata

type(t_mydata), target:: A(nSeries)

type(t_mydata), pointer:: Head

do i=1,nSeries-1

 A(i)%next => A(i+1)

end do

A(nSeries)%next => A(1)

Head=>A(1)

do while(something interrupts)

 call ReadData(Head%data)

 Head=>Head%next

end do

I doubt there would be performance difference between those two, but if you need to involve ALLOCATABLE arrays of A or A%data, you will change the memory layout, for the good or the bad of it: it might be more memory-fragmentation friendly or less cache-friendly. In any case, you should do your own benchmarking.

[jimdempseyatthecove]

Clemens,

Jugoslav had some good points.

Either a preallocatedarray or a linked listof arrays will work. There are some issues though to resolve regarding a "warmup" time between first sample and when the list of arrays become fully populated. Something along the lines of the following may be helpful

module foo 
 integer :: NumberOfArrays = 10000
 integer :: ArraySize = 10000

 integer :: Newest = 0
 integer :: Oldest = 0
 integer :: Insertion = 1
 real, allocatable :: Array(:,:)
 ...
 logical :: Done = .false.
end module foo

Program fooMain
 use foo
 ...
 allocate(Array(ArraySize, NumberOfArrays))
 call StartSampleThread
 do while(Oldest .eq. 0)
 sleep(100)
 end do
 do while(.not. Done)
 ! process valid arrays
 i = Oldest
 do while(i .ne. Insertion)
 call process(Array(:,i))
      if(i .eq. Insertion) call BufferOverflowed

 if(i .eq. Newest) exit ! finished
 i = mod(i + 1, NumberOfArrays) + 1
 end do
 end do
end Program fooMain

subroutine process(workinArray)
 use foo
  real :: workingArray(ArraySize) 
 ...
 do j=1,ArraySize
 ...
 end do
end subroutine process

subroutine InsertData
  use foo
 ... ! insert data
 Newest = Isertion
 if(Oldest .eq. 0) Oldest = Newest

 Insertion = mod(Insertion + 1, NumberOfArrays) + 1
 if(Insertion .eq. Oldest) Oldest = mod(Insertion + 1, NumberOfArrays) + 1
end subroutine InsertData

Jim Dempsey

[davinci]

Thanks guys.

(always these fast solutions from Jugoslav, Jim,Steve(Steveprobably with a day off on the beach)). I'll implement these methods and reply with the benchmark results asap.

Clemens (davinci)

[anthonyrichards]

Have you chosen your soubriquet so that everything you program is 'Davinci code' ?