Referencing external variables vs dummy variables within @einsum

[StephenVavasis]

I am in need of a package very much like Einsum.jl. Please refer to my julia-users posting, which is here:

https://groups.google.com/forum/#!topic/julia-users/CYGxXh1jNDg

as well as documentation of my vaporware, which is here:

https://github.com/StephenVavasis/indicial.jl

However, Einsum is missing a key feature that I need, namely, indirect addressing. In other words, I need to be able to transform a loop like

     x[:] = y[:] + z[i[:]]

into a for-loop.

Is there any possibility of adding indirect addressing?

Two other issues:

Broadcasting seems to work

   A[a,b] = x[a]  #OK

except in the case of no indices on the RHS, e.g.,

    x[a] = 0.0 #not OK

Also, the package does not seem to run under 0.5?

Thanks, Steve Vavasis

[ahwillia]

Hi Steve,

I'm spread a bit thin this week, but will try to address this soon. I'm not sure it would be an easy fix, but it would be nice to have this flexibility. Are there any use cases for having more than one level of indirect addressing? I.e. could we get away with just supporting z[i[:]] and not z[i[j[:]]]? Any other details you can provide on the applications of this would be useful.

[StephenVavasis]

It appears that I need subscripts nested only one deep. I am still typing in my code, and so far I have been writing it in terms of my nonexistent indicial package; here is an actual example of a statement with depth-1 nesting:

        @indicial xy[~1,~2] = bulkmesh.nodepos[mynodes[~1],~2]

I just found another case that doesn't seem to work in einsum. For this case, I want to make concatenated copies of a given integer vector:

        @indicial hesstmp[~1] = gradrows[mod(~1-1,ndof_elt) + 1]

I tried

    macroexpand(:(@einsum hesstmp[a] = gradrows[mod(a-1,ndof_elt)+1]))

but got an error message.

Thanks, Steve

[StephenVavasis]

If you can't get to this soon, I could still write my indicial package; it seems like einsum.jl would be a good starting point for indicial.jl since you've already figured out all the tricky issues. Would you mind if I copy your code? I'll give you credit in the comments & documentation.

[ahwillia]

It would be ideal if you could open a pull request to improve the code here. I'd greatly appreciate any improvements you make, as I will likely try to implement them later anyways.

If this is absolutely not possible, credit in the docs and comments is okay.

[StephenVavasis]

Alex,

It seems like the interface I need is incompatible with the one provided by einsum right now. In particular, I don't how in general the macro can distinguish the dummy indices a,b,c on the RHS from program variables. The proposal for @indicial solves this problem by using ~1, ~2, etc. for dummy indices instead of variable names. If you don't mind having two incompatible interfaces into the einsum package, then yes I'm happy to handle this with a PR.

[ahwillia]

We can use let blocks so that we don't overwrite index variables like i, j, and k. I just pushed a rough commit to fix this, see https://github.com/ahwillia/Einsum.jl/commit/d4202b111bcf069d11721ad16be40aacb07511f7

Now this new test passes:

i = -1
y = randn(10)
@einsum x[i] := y[i] 
@test i == -1 # rather than i == 10

The basic behavior of a let block is:

i = -1
let
    local i = 1 # var i in main namespace *not* overwritten
end
@show i  # prints i = -1

More documentation here: http://docs.julialang.org/en/release-0.4/manual/variables-and-scoping/

Thanks for the comments on the package. In general, I think we have pretty similar goals and I'd like to work synergistically as much as possible.

[StephenVavasis]

Alex,

I didn't make my previous message sufficiently clear. For my code, I need to support statements like this (simple example): x[~1] = y[~1,j] where j is a program variable and ~1 is a dummy index that should generate a for-loop. I don't see how to distinguish program variables from dummy indices on the RHS if the dummy indices look the same as variables. Possibly the macro could query the symbol table to determine that j has already been declared? I don't know how to do that.

[ahwillia]

There is probably a hacky way to do it. But I potentially like your solution better.

My philosophy for this package was that the index vars shouldn't be set like this. One way of getting around it would be to use subarrays:

Y = randn(10,10)
j = 3
y = sub(Y,:,j)
@einsum x[i] = y[i]

If this isn't satisfactory for you, then we could consider either having separate packages, or having two separate macros/interfaces in this repo. Would @indicial be effectively the same as @einsum other than this ~1, ~2 notation?

[StephenVavasis]

Alex,

Unfortunately, I cannot use subarrays for my purpose because the return value from the 'sub' function is created on the heap, not the stack. This means that a call to 'sub' in order to extract a small submatrix or subarray inside an inner loop is a performance-killer.

There are certainly advantages to using variable-name dummy indices rather than integers. One obvious advantage is that variable-name dummies can be made descriptive for better self-documentation. I'll try to think if there is some elegant way to combine the two interfaces.

[ahwillia]

What about @einsum x[i] = y[i,j] creates two dummy variables i and j, and sums over the second one. Whereas @einsum x[i] = y[i,:j] creates only one dummy variable i, and references a pre-defined variable j?

[StephenVavasis]

Alex,

I think that would work. For my example of making copies of a vector, I could say:

@einsum hesstmp[i] = gradrows[mod(i-1, :ndof_elt) + 1]

I wouldn't write :mod in place of mod because it is immediately followed by a parenthesis, so that macro can tell that it is a function name, right?

[ahwillia]

Yes, in principle I think that would be fine. Though we would need to add some code to get this working. I'd be interested in supporting this if you are still interested in opening a PR.

[StephenVavasis]

Alex, I created a branch, but for some reason I can't make a PR. I have a lot of trouble with git. Maybe you can turn my branch into a PR? The repository name is StephenVavasis/Einsum.jl and the branch name ref_external.

I implemented all of the features I need. I didn't make benchmarks or test cases yet, because I want to see if works in my real application. Unfortunately, I had to delete the capability to create the LHS using the := syntax. There were too many incompatibilities between my additions and the creation capability e.g., what should the following do?

  @einsum x[i,i] := y[i]
  @einsum x[i,:a] := y[i]

But I added += and -= capability.

[ahwillia]

Hmm did you create your version of Einsum.jl exactly? It looks like you didn't fork the repo? Maybe you created your github repo manually?

If you click the fork button on the upper right corner of this page, and then push your commits to that forked repo it should be straightforward to open a PR. I think the problem at the moment is that Github has no way of telling your repo is linked to this one (it thinks you made it from scratch). After you fork the repo you'll need to:

git remote add myremote https://github.com/StephenVavasis/Einsum.jl.git
git push myremote mybranch

In principle I don't think there should be a problem with using := but feel free to break things to get what you need working, we can fix/address that later.

[StephenVavasis]

Alex,

Here is how I created this repo: I started with a shell command 'git clone' to get a copy of your repo on my desktop. Then I made a branch on my desktop. I edited my branch, made a commit, and then I tried git push origin. This failed because I don't have write-access to your repo. So then I made a repo on my github account, and I carried out a git push to my repo on GitHub. However, there doesn't seem to be a useful button for creating a pull request back to your repo.

Anyway, I just followed your instructions above, and this seems to have fixed the problem, thanks!

[ahwillia]

Yeah the problem is that you can't make a new repo on your account, you need to fork it so that github knows what to compare it to.

Seeing as the conversation has diverged, I'm going to close this issue, and we can continue chatting on the pull request you opened: https://github.com/ahwillia/Einsum.jl/pull/3