EDIT (4/21/20): new link re: function inserted at end

Source: my tweets, with minor errors removed (https://twitter.com/krzhang/status/1252529588049072128)

Let's assume $A$ is symmetric and 0 on the diagonal. (disclaimer: I'm guessing this is not what you meant by "symmetry" because of poset context, but it may still be helpful) This means we are really working with unlooped undirected graphs, where
- p is how many edges the graph has, and
- q is now many 2-paths (disregarding order) that are not loops.

Now, cool observation: 2-paths that are not loops can be identified with their middle point and 2 neighbors. So for each vertex $i$ of deg. $d_i$, it contributes $d_i(d_i-1)/2$ 2-paths that are not loops.

So our problem becomes: "How many ways are there to split $p$ into nonnegative integers $d_1 + ... + d_n,$ such that $\sum d_i(d_i - 1)/2 = q$?"

Some manipulation gives $\sum d_i^2 = 2q + p$, so this problem really reduces to
"Given the first and 2nd power sums of $d_1 ... d_n$, how many sets of nonnegative $d_n$ are there?" or the quite-beautiful probabilistic form:

"How many nonnegative integral distributions are there of a fixed mean and variance?"

There're number theory constraints here, so I guess this is hard (which means original problem is even harder). However, computationally this isn't bad. Here's a solution:

- construct a 3-d infinite array so that $P[x][y][z]$ = "the number of ways to solve this with $x$ numbers such that their 1-power sum (sum) is $y$ and their 2-power sum is $z$
- use
*dynamic programming* to build this layer by layer by $x$. So compute everything with $x = 1$ first, then reduce each problem with $x+1$ to those with $x$ by summing over different values for the first element.

This gives an $O(n^2p(2q+p))$ algorithm.

Link (h/t Boris Alexeev): This last function is explored at https://mathworld.wolfram.com/SumofSquaresFunction.html as "sum of squares function." As I predicted, it seems number-theoretical and thus needs to be written as sums of modular functions for larger $n$ and $k$. Therefore, finding a closed form seems super hard.