Re: [ferret_users] count consecutive events

[mick spillane <mick.spillane@xxxxxxxx>]

Hi Patrick,

   I had a similar problem some years ago - counting runs in sequences 
of 0's and 1's -- but just thought of a better way to do it, using the 
@cia transformation.  One identifies "start" and "end" points in the 
sequence then, whereever a "start" is found, looking for where the "next 
end above" is located.

   The method results in short, independent definitions for each run 
length, and lends itself to use in a repeat loop.  A demo script (when 
the sequence is along the x-axis) is attached but the key points are:

! define start points (a 1 preceded by a 0) for the runs in the sequence "v"
let/title="Run Starts" start=if(v eq 1 and v[i=@shf:-1] eq 0)then 1
! ... and end points (a 1 followed by a 0)
let/title="Run Ends" end=if(v eq 1 and v[i=@shf:1] eq 0)then 1
! then the start points of runs of length "n" are given by the mask "run"
let run=IF(START EQ 1 AND END[I=@CIA] EQ `n-1`)then 1
! and to count them
list run[i=@ngd]

Hope it helps,
Mick
---------
Brockmann Patrick wrote:
> Hi all,
>
> Is there a way to count consecutive events ?
> From climate indice variables, I would like to have the number
> of continuous events, for 1, 2, 3, 4, ..., n events.
>
>               B   C1EVN  C2EVN  C3EVN  C4EVN
> 1    /  1:  0.000   ....   ....   ....   ....
> 2    /  2:  1.000  1.000   ....   ....   ....
> 3    /  3:  0.000   ....   ....   ....   ....
> 4    /  4:  1.000  1.000   ....   ....   ....
> 5    /  5:  0.000   ....   ....   ....   ....
> 6    /  6:  1.000  1.000   ....   ....   ....
> 7    /  7:  0.000   ....   ....   ....   ....
> 8    /  8:  1.000  1.000   ....   ....   ....
> 9    /  9:  0.000   ....   ....   ....   ....
> 10   / 10:  0.000   ....   ....   ....   ....
> 11   / 11:  0.000   ....   ....   ....   ....
> 12   / 12:  1.000   ....  1.000   ....   ....
> 13   / 13:  1.000   ....   ....   ....   ....
> 14   / 14:  0.000   ....   ....   ....   ....
> 15   / 15:  1.000   ....  1.000   ....   ....
> 16   / 16:  1.000   ....   ....   ....   ....
> 17   / 17:  0.000   ....   ....   ....   ....
> 18   / 18:  0.000   ....   ....   ....   ....
> 19   / 19:  1.000   ....  1.000   ....   ....
> 20   / 20:  1.000   ....   ....   ....   ....
> 21   / 21:  0.000   ....   ....   ....   ....
> 22   / 22:  0.000   ....   ....   ....   ....
> 23   / 23:  1.000  1.000   ....   ....   ....
> 24   / 24:  0.000   ....   ....   ....   ....
> 25   / 25:  0.000   ....   ....   ....   ....
> 26   / 26:  1.000   ....   ....   ....  1.000
> 27   / 27:  1.000   ....   ....   ....   ....
> 28   / 28:  1.000   ....   ....   ....   ....
> 29   / 29:  1.000   ....   ....   ....   ....
> 30   / 30:  0.000   ....   ....   ....   ....
> ...
>
> C1EVN represents events for 1 period at 1
> C2EVN represents events for 2 consecutive periods at 1
> C3EVN represents events for 3 consecutive periods at 1
> C4EVN represents events for 4 consecutive periods at 1
> ...
>
> Attached is a exemple I have scripted. It works well for a 52 points 
> variable
> when I look for for 1 to 4 consecutive events.
>
> As the real case is to compute this over a year of 360 days looking for
> consecutive events of 30 days, I need to define variables from very 
> long expressions.
> and in this case I get the message :
> **ERROR: expression too complex
>
> Let me know if there a solution to this.
> Patrick
åå

! Counting runs of 1's in a sequence of 0's and 1's

! generate a series of 0's and 1's and list to a file ...
let rv=if(randu(i) ge 0.45)then 1 else 0
list/nohead/i=1:360/file=runtest.d/form=(f3.0) rv
! ... then read them back
def axis/x=1:360:1 xax ; def grid/x=xax grd
file/form=free/var=v/g=grd runtest.d

! define start points for the runs (a 1 preceded by a 0)
let/title="Run Starts" start=if(v eq 1 and v[i=@shf:-1] eq 0)then 1
! ... and end points (a 1 followed by a 0)
let/title="Run Ends" end=if(v eq 1 and v[i=@shf:1] eq 0)then 1

! define some simple runs (of 1, 2, 3, or 4 ones in a row)
! a run of length 1 starts and ends at the same point
! let run1=if(start*end eq 1)then 1
! ... or equivalently
let run1=if(start eq 1 and end[i=@cia] eq 0)then 1
! a run of length 2 has its end "one point ahead" ...
let run2=if(start eq 1 and end[i=@cia] eq 1)then 1
! a run of length 3 has its end "two points ahead" ...
let run3=if(start eq 1 and end[i=@cia] eq 2)then 1
! and so on
let run4=if(start eq 1 and end[i=@cia] eq 3)then 1
! since these expressions are short and independent you can
! define as many as you like ...

! ... and verify that you match the results from a test listing
list/nohead/i=1:52/form=(f2.0) v
list/nohead/i=1:52 run1[i=@ngd],run2[i=@ngd],run3[i=@ngd],run4[i=@ngd]

message

! to look at multiple run lengths you could employ a repeat loop
! --- lets increase the probability of "ones" to get some longer runs

let rv=if(randu(i) ge 0.1)then 1 else 0
list/nohead/i=1:10000/file=runtest.d/form=(f3.0)/clob rv
can data 1 ; def axis/x=1:10000:1 xax ; def grid/x=xax grd
file/form=free/var=v/g=grd runtest.d

message

! lets list the number of runs of lengths 1-30 in this 10000 point sequence
sp touch run.results
repeat/range=1:30/name=n (let run=IF(START EQ 1 AND END[I=@CIA] EQ `n-1`)then 1 ;\
list/nohead/app/form=(2f6.0)/file=run.results n,run[i=@ngd] )

! and make a bar chart of the results
let bx=ysequence({-0.5,0.5,0.5,-0.5}) ; let by=ysequence({0,0,1,1})
file/form=free/var=length,hits run.results
poly/nolab/pal=black/coord=y/set bx+length,by*hits
ppl yfor,(i3) ; ppl xfor,(i2) ; ppl xaxis,0.5,30.5,1 ; ppl axlint,1 ; ppl fillpoly

q