Generating variables using indexing

[Jason Grout]

This question has come up often enough that I thought posting a short
snippet here might be useful. The question is: how do I automatically
generate variables based on indices (e.g., a[0], a[1], etc. being
variables). Here is one way:

class VariableGenerator(object):
def __init__(self, prefix):
self.__prefix = prefix

@cached_method
def __getitem__(self, key):
return SR.var("%s%s"%(self.__prefix,key))

Now just specify a prefix, and then you can index to your heart's
content to generate variables.

a=VariableGenerator('a') # some people may like 'a_' as the prefix
a[0], a[1], a[2] # all variables

Of course, this can easily be extended to using function call syntax:
a(0), or to using multiple indices: a[1,3]. Indeed, you can let your
imagination run wild and even do things like return full symbolic
matrices or vectors with slices: a[0:5, 0:5].

Perhaps this is useful enough to be in Sage instead of just a snippet
here too...

Thanks,

Jason

--
Jason Grout

[kcrisman]

Uh, yeah.

With the proviso that the usual answer to the question should also be
available, so that if people want

a1, a2, a3, a4

instead of

a[1],a[2],a[3],a[4]

that is an option too. I was wondering about how hard it would be to
make the former easily available a few weeks ago, but didn't follow
through.

See http://trac.sagemath.org/sage_trac/ticket/11576 for adding this.

[GaryMak]

Hi guys

apologies as usual for how dumb this question's going to sound ... but how do I actually use the variables a[1] etc? What I was hoping to do was to use them in the following way (once I had already apparently successfully invoked Jason's routine), where d is a variable I would like to be able to change:

d = 5;

for ii in range(0,d):

    for a[ii] in range(0,d):

        print a[ii]^2 + ii;

but I get "TypeError: 'VariableGenerator' object does not support item assignment" seemingly no matter what variant I try.

Thanks for any pointers!

kind regards

Gary

[Tom Boothby]

Gary,

The third line

for a[ii] in range(0,d):

should read

for ii in range(0,d):

> --
> You received this message because you are subscribed to the Google Groups
> "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-support...@googlegroups.com.
> To post to this group, send email to sage-s...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-support?hl=en.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Tom Boothby]

oops, didn't read carefully enough... that third line just doesn't
make much sense.

[Gary McConnell]

I guess that's what I don't understand then: why can I not invoke a structure like this? In fact I need to nest several levels of this type of structure, so I would then define variables b_xx etc as b[a[ii]] and so on. If that is not how these variable definitions work, could someone please tell me how they do work? ie in what situation would I be able to use the a[0], a[1] etc defined in Jason's email, as variables (either as loop counters or for anything else) - what syntax is required?

thanks

You received this message because you are subscribed to a topic in the Google Groups "sage-support" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sage-support/GFJdjFvKCvo/unsubscribe?hl=en.
To unsubscribe from this group and all its topics, send an email to sage-support...@googlegroups.com.

[Jason Grout]

On 3/29/13 6:38 PM, GaryMak wrote:
> d = 5;
> for ii in range(0,d):
> for a[ii] in range(0,d):
> print a[ii]^2 + ii;
>

What do you expect the output of this to be? What are you trying to
accomplish with this code?

Thanks,

Jason

[Gary McConnell]

Hi Jason

fair question - I was trying to avoid including all the boring details of my code! In the simplest case I have a search routine which looks for M sets of d vectors inside a vector space of dimension N over a finite field, whose dot products satisfy a bunch of polynomial equations. My problem is that I need to compare each new vector with all of the antecedent vectors, so that the number of "equations" varies with d, M and N (in several "recursive dimensions"). The only way I have been able to do it so far is to construct separate loops by hand for each of d,N,M and for each "new" vector within them ... ie every time I alter any one of d, N or M, I need a new program!!

So I had hoped that your code would have allowed me to "telescope" the search and comparison loops as functions of d,M,N in a neat way; only obviously I've misunderstood what it does. The ii, jj, kk, "a[ii,jj,kk]" etc are ONLY used as indexing variables in my hoped-for program; they are never invoked as anything else. I then have vectors "vex[a[ii,jj,kk]]" etc which I compare with one another etc etc. Note I have put the indexing in inverted commas because I hadn't figured out yet how to do >1 level of recursion inside your code still ....

I hope that clarifies somewhat what I am trying to do - thanks a lot

Gary

--
You received this message because you are subscribed to a topic in the Google Groups "sage-support" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sage-support/GFJdjFvKCvo/unsubscribe?hl=en.

To unsubscribe from this group and all its topics, send an email to sage-support+unsubscribe@googlegroups.com.

[Jason Grout]

On 3/30/13 5:08 AM, Gary McConnell wrote:
> Hi Jason
>
> fair question - I was trying to avoid including all the boring details
> of my code! In the simplest case I have a search routine which looks for
> M sets of d vectors inside a vector space of dimension N over a finite
> field, whose dot products satisfy a bunch of polynomial equations. My
> problem is that I need to compare each new vector with all of the
> antecedent vectors, so that the number of "equations" varies with d, M
> and N (in several "recursive dimensions"). The only way I have been able
> to do it so far is to construct separate loops by hand for each of d,N,M
> and for each "new" vector within them ... ie every time I alter any one
> of d, N or M, I need a new program!!
>
> So I had hoped that your code would have allowed me to "telescope" the
> search and comparison loops as functions of d,M,N in a neat way; only
> obviously I've misunderstood what it does.

I think the confusion may be in what we mean by "variable". My code in
this thread is for easily constructing Sage symbolic variables, and is
just a convenient way to do something like var('x1,x2,x3') (x[1] is
something like var('x1')). It sounds like you are referring to python
variables, which are a different concept.

Still, if you can give a short example illustrating exactly what you
trying to do, I think we can still probably point you in the right
direction (for example, give us a working example of the nested for
loops that you currently have that you'd like to generalize). If I
followed you correctly, you should be able to do what you want fairly
easily.

Thanks,

Jason

[Gary McConnell]

Hi Jason - thanks a lot - I have attached a text file which is just a cleaned-up version of my code, with some explanations at the top for functions whose details are not relevant to the question. Apologies - you'll see it's a "neat mess", so to speak ....

Kind regards

Gary

Thanks,

Jason

[Jason Grout]

On 3/30/13 6:34 AM, Gary McConnell wrote:
> Hi Jason - thanks a lot - I have attached a text file which is just a
> cleaned-up version of my code, with some explanations at the top for
> functions whose details are not relevant to the question. Apologies -
> you'll see it's a "neat mess", so to speak ....
>

So to summarize your code, you're finding M=3 sets of d=3 mutually
orthogonal vectors so that each pair of sets satisfies some criteria
(your XRY function)? Do the sets have to be mutually exclusive? Your
code doesn't enforce this, but your xRy function may imply it. Also, if
you have an isotropic vector v (orthogonal to itself), then I suppose
your code allows the set [v,v,v]. Is that what you want?

Here's an example of how to rethink at least part of the problem:
finding a set of d distinct mutually orthogonal vectors:

N = set(N1) # copy N so we can modify it
for v in N:
base1 = [v]
N.remove(v) # we already have v
for w in N:
if all(w*vv==0 for vv in set1):
base1.append(w)
N.remove(w)
if len(base1)==d:
break
else:
# we went all the way through N, but didn't find enough vectors
raise ValueError("couldn't find enough mutually orthogonal
vectors")

Thanks,

Jason

[Jason Grout]

On 3/30/13 7:12 AM, Jason Grout wrote:
> So to summarize your code, you're finding M=3 sets of d=3 mutually
> orthogonal vectors so that each pair of sets satisfies some criteria
> (your XRY function)? Do the sets have to be mutually exclusive? Your
> code doesn't enforce this, but your xRy function may imply it. Also, if
> you have an isotropic vector v (orthogonal to itself), then I suppose
> your code allows the set [v,v,v]. Is that what you want?

Here's a more general function that I think (I haven't run the function)
finds *a* list of sets that I think satisfy your criteria. I think your
code finds all such sets, though. Is that what is really needed?

https://gist.github.com/jasongrout/5276508

Thanks,

Jason

[Gary McConnell]

Hi

in fact the construction of N1 precludes isotropic vectors - I should have mentioned that. The xRy function effectively forces separate sets of vectors to be mutually exclusive also. I have learnt a lot from your example there - but unfortunately I cannot get to where I need to - I have just received your new code ... will see how close we get - you are right in that I am only looking for one example - in fact the whole code is wrapped in the break function thing you showed me at https://groups.google.com/forum/?hl=en&fromgroups=#!searchin/sage-support/garymakonel/sage-support/IKChSZ9zdsY/rqeXHjusPxsJ !! and so I exit upon finding one complete set ...

Thanks,

Jason

[Gary McConnell]

Hi again - thanks for that excellent new function - it's amazingly compact and I think it may well do what I need - I am about to try modifying it to work in my context, and will let you know.

Am I to deduce from this that what I was originally trying to do (ie telescoping recursive variables) is not possible in Python?

Thanks,

Jason

[Jason Grout]

On 3/30/13 7:49 AM, Gary McConnell wrote:
> in fact the construction of N1 precludes isotropic vectors - I should
> have mentioned that. The xRy function effectively forces separate sets
> of vectors to be mutually exclusive also. I have learnt a lot from your
> example there - but unfortunately I cannot get to where I need to - I
> have just received your new code ... will see how close we get - you are
> right in that I am only looking for one example - in fact the whole code
> is wrapped in the break function thing you showed me at
> https://groups.google.com/forum/?hl=en&fromgroups=#!searchin/sage-support/garymakonel/sage-support/IKChSZ9zdsY/rqeXHjusPxsJ

> <https://groups.google.com/forum/?hl=en&fromgroups=#%21searchin/sage-support/garymakonel/sage-support/IKChSZ9zdsY/rqeXHjusPxsJ>

> !! and so I exit upon finding one complete set ...

I just noticed a problem with my code, and put a comment on the gist:
https://gist.github.com/jasongrout/5276508

Jason

[Jason Grout]

On 3/30/13 7:58 AM, Gary McConnell wrote:
> Am I to deduce from this that what I was originally trying to do (ie
> telescoping recursive variables) is not possible in Python?

Sorry---it wasn't ever really clear to me how you were approaching the
problem and what you meant by "telescoping recursive variables".

Thanks,

Jason

[Gary McConnell]

yes indeed that is a problem in the general case but in fact in my experience with the particular things I am looking at is that either the first such set works, or else none of them do .... so it may not be as much of a problem in this case. Fingers crossed!

Jason

[Gary McConnell]

fair enough - I'm not too clear myself on how to describe it sorry -  but for example, if you look at the code I posted in the .txt file, ideally the vv_i,ww_i,xx_i,yy_i etc would be variables which would only be created if the requisite antecedent vectors had been found (as effectively happens in my code) - then suddenly they would "telescope" out into a new range of d variables for the next level of the search.

Thanks,

Jason

[Gary McConnell]

One more thing Jason: the code you wrote, as far as I can see, does not take account of "N", which I guess is your N=set(N1) of my "N1" - the point is I guess that if I pass N1 to the function instead of your F then it should be ok (ie it only searches within the subset N1 of the full vector space); but (and this is how I got into this mess in the first place!) I keep on getting messages about mutable vectors being unhashable, no matter how I try to weave the list N1 into the search function. Hence my klunky kode ...

So perhaps I have inadvertently created N1 in such a way that it will not let me copy it; or is there another way to get that type of assignment "N=set(N1)" to work?

[Jason Grout]

On 3/30/13 8:39 AM, Gary McConnell wrote:

> One more thing Jason: the code you wrote, as far as I can see, does not
> take account of "N", which I guess is your N=set(N1) of my "N1" - the
> point is I guess that if I pass N1 to the function instead of your F
> then it should be ok (ie it only searches within the subset N1 of the
> full vector space); but (and this is how I got into this mess in the
> first place!) I keep on getting messages about mutable vectors being
> unhashable, no matter how I try to weave the list N1 into the search
> function. Hence my klunky kode ...
>
> So perhaps I have inadvertently created N1 in such a way that it will
> not let me copy it; or is there another way to get that type of
> assignment "N=set(N1)" to work?

Yes, you can pass in N1 instead of F --- all it does is iterate through
that set.

I added this line to make the vectors immutable before storing them in
the set: https://gist.github.com/jasongrout/5276508#file-findsets-py-L22

Thanks,

Jason

[Gary McConnell]

Brilliant - thanks Jason - sorry for delay

Jason

[Gary McConnell]

OK now I'm back in a familiar nightmare - if I set the vectors to be immutable then I cannot treat them as vectors because they seem to be of "NoneType" or something so have no length etc etc; however I cannot not set them to be immutable if I wish to work through a list consisting of those very vectors. I remember now that this is what drove me to abandon anything other than the plain old nested subroutines a long time ago, because I couldn't reverse-engineer my way out of this conundrum. Is there something simple that will put me out of my misery please?!!

[Jason Grout]

On 3/30/13 1:43 PM, Gary McConnell wrote:
> OK now I'm back in a familiar nightmare - if I set the vectors to be
> immutable then I cannot treat them as vectors because they seem to be of

> "NoneType" or something so have no length etc etc; however I cannot /not
> /set them to be immutable if I wish to work through a list consisting of

> those very vectors. I remember now that this is what drove me to abandon
> anything other than the plain old nested subroutines a long time ago,
> because I couldn't reverse-engineer my way out of this conundrum. Is
> there something simple that will put me out of my misery please?!!
>

Sure; convert to tuples to store into the set and test membership in the
set. I've updated the gist:

https://gist.github.com/jasongrout/5276508

Thanks,

Jason

[Gary McConnell]

Thanks again Jason - that seems to do the trick - now I just have to figure out the equivalent of the old BASIC "goto 10" to start your loop again!!

So I understand something of what is going on, could you please tell me why the tuple thing is necessary for the "used" set but the vectors are necessary for the "bases": ie why can we not use a vector structure for "used"?

Thanks,

Jason

[Gary McConnell]

OK so thinking about it, even though your code is beautifully compact and elegant, I think I am going to have to revert a little to the "outer-inner-loop" structure in order to achieve what I need. Namely, instead of storing "used" vectors, I store "used" bases and search through the remaining orthogonal sets exhaustively. Otherwise I cannot get at the bases which may share a vector with another basis but have different "xRy" properties. I am trying to program this now - please feel free to tell me a better way!!

Thanks and regards

Gary

[Jason Grout]

On 3/30/13 4:09 PM, Gary McConnell wrote:
> So I understand something of what is going on, could you please tell me
> why the tuple thing is necessary for the "used" set but the vectors are
> necessary for the "bases": ie why can we not use a vector structure for
> "used"?

The bases is just a list of lists of vectors, which doesn't require
immutability. The used variable is a set (since set membership testing
is O(1), but list membership testing is O(n)---a lot slower). A set
requires immutability.

Does that help?

Thanks,

Jason

[Jason Grout]

On 3/30/13 4:56 PM, Gary McConnell wrote:
> OK so thinking about it, even though your code is beautifully compact
> and elegant, I think I am going to have to revert a little to the
> "outer-inner-loop" structure in order to achieve what I need. Namely,
> instead of storing "used" vectors, I store "used" bases and search
> through the remaining orthogonal sets exhaustively. Otherwise I cannot
> get at the bases which may share a vector with another basis but have
> different "xRy" properties. I am trying to program this now - please
> feel free to tell me a better way!!

Not a problem. I was hoping my code was a good starting off point to
rethinking the problem. I was playing with a way to just have a
function that generated all possible orthogonal sets of d vectors each,
and then you could test your relationship for all of them in all
possible ways. Maybe that's the best thing.

However, maybe at this point it would be good to ask on, say,
sage-combinat mailing list. You have a very combinatorial problem here,
and they have lots of tools for enumerating sets with various
constraints. For example, looking at
http://www.sagemath.org/doc/reference/combinat/sage/combinat/backtrack.html,
you could use SearchForest to implement your "find all independent sets
of size d" pretty easily. The roots would be the list of singletons,
and the children of a list could be found by:

* if the list has d elements in it, it has no children
* if the list has less than k<d elements in it, then you return all
lists of size k+1 where you basically append all vectors that are
orthogonal to all the things in your list.
* to make this efficient, you want to iterate through your vectors in
some sort of order, and then only test and add vectors that are ordered
past your last vector in your current list.

For example:

V=FiniteField(3)^4
L = list(V)
d=9
def find_children(node):
if len(node)==d:
return []
node_position = L.index(node[-1])
if len(node)+len(L)-node_position < d:
# not enough elements to finish off the list
return []
children = []
# find the last element
for v in L[L.index(node[-1])+1:]:
if all(v*w==0 for w in node):
children.append(node+[v])
return children

def post_process(node):
# only output nodes with the right length
if len(node)==d:
return node
else:
return None

s=SearchForest(roots=[[v] for v in L], children=find_children,
post_process = post_process,
category=FiniteEnumeratedSets())

# mutually orthogonal sets of size 9
A=list(s.depth_first_search_iterator())
print len(A)
print A

(see http://sagecell.sagemath.org/?q=e2a9c6b7-bb0c-40e6-9081-3609aaa3227f)

Anyway, s.depth_first_search_iterator() gives an iterator over all
independent sets of size d. Now you can test for your relations between
the independent sets. In fact, you could set up another SearchForest
for that too :). The nice thing about SearchForest is that it takes
care of all the pruning and other things for you. You can do
breadth-first or depth-first searches in it, etc.

Another way to think about your problem is: you have a graph with
vertices = vectors in N1, where two vectors are connected if their dot
product is nonzero. You're trying to find a list of M mutually
exclusive independent sets of size d, where the independent sets must
also satisfy some relationship relative to each other. Of course, you
could also take the complement and look for cliques. The combinat
people might have already implemented fast searching for independent
sets, etc.

Feel free to share your code with us; we might be able to give you
further pointers. It'd also be fun to see where this research is headed
(e.g., a link to a paper), if it's appropriate. It'd be cool to see
another paper in the Sage library links:
http://www.sagemath.org/library-publications.html

Thanks,

Jason

[Gary McConnell]

Hi again

well, not that I understand how, but that code works magnificently! I see (confusingly for me I am afraid :) ) that the concept of "reserved words" does not have much currency in SAGE ... ! 

The only difficulty with using that code is that for the sets I am looking at, the number of bases in the tree grows horribly quickly with p,d etc and so I still will need a routine which "checks as I go along", rather than one which stores everything first and then searches. I tried and failed to write a forest thing for that, then wrote a dumb nested loop thing which takes as a starting point your tree of nodes, and verified that in "small" cases I do indeed get the right results overall. But it is horrendously slow, presumably because of the memory requirements for storing large-ish sets of bases.

So somehow I need to understand how to incorporate the XRY condition (ie X and Y as sets of vectors are pairwise related by the R condition - ie xRy for all x in X and for all y in Y) as I crawl along the tree - then we stop if we find M examples (ie a node with level M). Presumably knowing "not M" will be pseudo-equivalent to the Halting Problem, sort-of, so I am not concerned about that!

Jason I have sent you a separate email regarding the points you made in your last post - I just need a spot of guidance there - but briefly the xRy condition is as follows: my original search space consists (essentially entirely with some boring caveats) of the subSET of vectors of F^N whose "norm" (ie dot product with itself) is 1 and all of whose entries come from a fixed (finite!) subset S_d = S(d,F) of the underlying field F, which in particular does not contain 0 or 1.

Observe by the way that d<=N by definition (I hadn't thought of relaxing the norm<>0 condition and taking sets with d>N as you did - perhaps some new research lies there!!). Then xRy iff x.y \in S_d, where x.y can for now be taken to be the dot product. Also if v \in X and v \in Y then XRY CANNOT HOLD (since v.v=1 and 1 is not in S_d). So any bases which have a vector in common do not need to be tested; similarly any pairs of bases containing one from each which are orthogonal to one another can be discarded (ie v \in X, w \in Y and v.w=0 but 0 not in S_d so NOT XRY).

So one needs to check that the "Gram matrix" of each pair of bases {X,Y} contains only entries from S_d; if so, then XRY. My "code" consists of nothing more than setting up S_d and checking these conditions in a very obvious way. 

Thanks again for the marvellous help so far.

Best regards

Gary

Thanks,

Jason

[Gary McConnell]

One other thought - if we do manage to get this search routine to work "sequentially" (ie not necessarily to have to hold all things in memory at the same time), then perhaps we can extend this philosophy right back to the definition of "N_1" (my candidate vectors with norm 1 etc), because once again for p>5 and d>3 and N>=d even the size of N_1 becomes so large that Jason's algorithm also grinds to an effective halt in subsequently trying to generate all of the orthogonal "bases" (they are bases obviously iff d=N, but let me refer to them all as bases for simplicity), let alone trying then to search through them ... I should mention that a simple rule of thumb for the size of our search set of vectors N_1 is that it is of order something like q^((N+1)/2) in low dimensions, where q is the size of the underlying field.

So I need to begin to construct and test N_1, most fundamentally, and while doing that "catch" each successive d-tuple of vectors forming a basis, and store it in its own parallel tree, and then while doing that in turn catch each successive 2-, 3-, 4- tuple of such bases {B,C,D,...} satisfying B_R_C and B_R_D and ... and C_R_D and ....  It is in these final tuples of bases of possibly varying length that we are ultimately interested; in my original posting I mentioned a "holy grail" number M = M(q,d,N) which is the maximum size of the last type of tuple of bases - eg it is a big deal (in some cases!) if we can find a set satisfying M=d+1.

But in essence the point I am trying to get across is that I do not care at all what N_1 looks like, nor do I care about the bases, other than if they constitute one of the growing sets of bases at the end of the "tree" - ie we may discard all "dead-end" information as we go along.

Would it perhaps help for me to post a simplified meta-code version of this algorithm? The unnecessary details buried in the actual code are very distracting ...

[Gary McConnell]

yes, thanks!

Thanks,

Jason