Number of acyclic digraphs (DAGs) with n labelled nodes and r arcs

Question

I am trying to find a way to calculate the number of acyclic digraphs with $n$ labelled nodes and $r$ arcs.

In his paper, Robinson defined a counting function for acyclic digraphs and on the final page of this paper, under 3. Additional Results, he says that "standard is the introduction of the total number of arcs as an additional parameter". Sadly, I'm a bit of an idiot, so introducing this additional parameter is a a bit beyond my understanding.

I have found that Rodionov here (page 320, following equation (4.) under Theorem 1) finds the following recurrence relation for number of DAGS, $A_{n,r}$, where $n$ is number of nodes and $r$ is number of arcs, as I would seem to like. His result appears to be as follows:

$A_{n,r} = \sum_{m=1}^{n} (-1)^{m-1} \binom{n}{m} \sum_{k=0}^{r} \binom{m(n-m)}{r-k} A_{n-m,k}$.

Is this result correct? I have tried calculating out a few values by hand to test it. I'm unsure if I am reading the formula wrong, if I am misunderstanding something fundamental or if my calculations are incorrect - but something is going amiss.

I used the following values to test the above recurrence relation: $A_{0,0} = 1$ $A_{1,0} = 1$ $A_{2,1} = 2$ $A_{2,0} = 1$ $A_{3,0} = 1$ $A_{3,1} = 6$ $A_{3,2} = 9$, which didn't seem to work.

Another issue, the recurrence relation means you have to define values like $A_{1,1}$, $A_{0,1}$ and $A_{0,0}$, where are nonsensical. How can you have an arc without a nodes? How can you have a node and an arc - unless you allow self-loops? Do we just set these terms to $0$ and forget about them?

If you could write a python function that spat out the correct answer for a given $n$ and $r$, you are a much cleverer person than me.

Edit:

I have found the OEIS(A081064) link here, which makes me feel confident that the formula is correct, since it appears to have been successfully implemented in Maple - although I do not understand the formula.

Please see my Python script below, since it was not included before, which I would still like to get to work:

from scipy.special import comb

def sum_1(n,r):
    output = 0
    for m in range(1,n):
        output = output + ((-1)**(m-1))*comb(n,m)*sum_2(n,r,m)
    return output

def sum_2(n,r,m):
    output = 0
    for k in range(0,r):
        output = output + comb((m*(n-m)),(r-k)) * DAGs(n-m,k)
    return output

def DAGs(n,r):
    if n==1 and r == 0:
        return 1
    else:
        return sum_1(n,r)

Also, please see a list of correct values for reference: $A_{1,0} = 1$ $A_{2,0} = 1$ $A_{2,1} = 2$ $A_{3,0} = 1$ $A_{3,1} = 6$ $A_{3,2} = 12$ $A_{3,3} = 6$ $A_{3,4} = 6$ $A_{4,0} = 1$ $A_{4,1} = 12$ These have been validated by hand and against the OEIS link.