Antecedent Precipitation Index
Vector 13:3, page 8
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 8
In his recent letter (Vector Vol.13 No.2, "While in J …") Nick Cox asks whether a quicker concise version can be produced of his J verb "api2", which is reproduced here:
api2 =. 3 : 0
NB. antecedent precipitation index NJC 15 April 1996
NB. k api p
:
+/"1 (* x. &^ @ i. @ - @ #) \ y.
)
0.5 api2 1 2 3 4
1 2.5 4.25 6.125
It can be seen that, given a vector p and a scalar k, then (k api2 p) computes the vector a from the relationship (which I restate as):
a[0] = p[0]
a[i+1] = p[i+1] + k × a[i]
Nick Cox notes that api2 is slower and needs more space than his longer version:
api1 =. 3 : 0
NB. antecedent precipitation index NJC 15 April 1996
NB. k api p
:
k =. x.
p =. y.
api =. 0
i =. 0
while. i < #y. do.
api =. api , (i { p) + k * {: api
i =. i + 1
end.
}. api
)
The longer version api1 implements the recurrence relation directly, while the shorter api2 implements the derived relation:
that is, each item of a is a polynomial in k. The two versions give the same results:
Performance of the two versions was compared, at 3 different lengths for the vector argument p (with p=. i.10, p=. i.100 and p=. i.1000). Times in seconds on my machine, at the 3 argument-lengths, are:
api2 gives "out of memory" at argument length 1000, and is slower at argument length 100. The first improvement which might be made is to do the summing as soon as possible, rather than leaving it to the end. This means that we can avoid building an n by n matrix as an intermediate result. The verb I differs from api2 only by absorbing the +/ into the verb which is given to the \ operator:
The timings for the three versions are:
Evidently I is no faster than api2, but it does run to completion. The next improvement can be made by using J’s primitive p. to evaluate all the polynomials, rather than the programmed equivalent of api2. The algorithm is essentially the same, although reversals of ordering are needed to conform with p. .
Verb II is faster than I, although still not as good as api1 at the longer argument lengths:
Even so, algorithms api2, I and II are inherently inefficient because each term in the result is computed independently: all that work is thrown away when with just a little more work it could be used to give the next term as well! In relation to the argument-length, the execution time of api2, I and II looks like a square-law, whereas that of api1 looks linear.
Let us return to the serial algorithm of api1. If we use a J operator to replace the programmed control of api1, the iterator ^: (power conjunction) or the insertion operator / both look promising. Both were tried, but / gave the better result.
With / we will need a dyadic verb, u say, to insert between each item of argument p. Given a value for k:
then a first approximation to u would be:
which closely resembles the key line of api1. Inserting u between each item of (the reverse of) argument p gives the required result:
Now we need to supply k as an argument rather than as a global variable. k can’t be an argument to the dyad u because both of u’s argument-places are spoken for. If we first define verb v as a tacit equivalent to u:
then we can parameterise v to make an adverb, A, say, and the argument of A will be the value k.
Combining A with the insertion operator, and the needed reversal of the argument p we get an adverb, III say, which is structured as a train of 3 adverbs:
The timings of all versions are:
10 100 1000
api1 0.05 0.43 4.35
api2 0.02 0.53 _
I 0.02 0.52 59.03
II 0.01 0.14 9.30
III 0.02 0.09 1.37
Algorithm III is faster than api1 at all argument lengths tried.