bapolar.texw

\documentclass{amsart}
\usepackage{amsmath,amssymb,amsthm,mathrsfs,microtype}
\usepackage{minted}
\usepackage[margin=1in]{geometry}

\newtheorem{theorem}{Theorem}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{lemma}[theorem]{Lemma}

\newcommand{\CC}{\mathbb{C}}
\newcommand{\ot}{\otimes}

\newcommand{\abs}[1]{\lvert #1 \rvert}
\newcommand{\Abs}[1]{\left\lvert #1 \right\rvert}
\newcommand{\norm}[1]{\lVert #1 \rVert}

\newcommand{\code}[1]{\mintinline{python}{#1}}

\DeclareMathOperator{\LT}{LT}
\DeclareMathOperator{\LM}{LM}
\DeclareMathOperator{\Lm}{\textsc{Lm}}
\DeclareMathOperator{\ltcoeff}{\textsc{LtCoeff}}
\DeclareMathOperator{\wt}{wt}

\title{Border apolarity implementation}
\author{Austin Conner}

\begin{document}
\maketitle
\section{Introduction}\label{sec:intro}
Given a tensor $T \in A \ot B\ot C$, a natural number $r$, and a solvable
connected Lie group $B \subset \operatorname{GL}(A) \times \operatorname{GL}(B)
\times \operatorname{GL}(C)$ stabilizing $T$, the routines of this file seek to
practically implement complete enumeration of multigraded ideals $I\subset
\operatorname{Sym}(A^* \oplus B^* \oplus C^*)$ satisfying
\begin{enumerate}
\item \label{Bfixed} I is $B$-fixed
\item \label{inTperp} $I \subset T^\perp$
\item \label{codim} The $(ijk)$-graded component of $I$, $I_{ijk} \subset S^i(A^*)\ot S^j(B^*)
\ot S^k(C^*)$, has codimension $\min (r, \dim S^i(A^*)\ot S^j(B^*) \ot
S^k(C^*))$.
\end{enumerate}

Such ideals may occur in positive dimensional families, so more precisely we
wish to enumerate computational descriptions of families of ideals which
together exhaust all ideals satisfying the conditions. The kinds of families of
ideals manipulated by the program are described in \S\ref{sec:representation}
and called \emph{representable}. A representable family of ideals $I_t$ always
satisfies conditions (\ref{Bfixed}) and (\ref{inTperp}) above, and $t$ ranges
over an affine variety. Hence, the work of the program is to enumerate an
exhaustive list of representable families satisfying condition (\ref{codim}). 
The program proceeds by exhaustively enlarging ideals so that condition (\ref{codim}) is
satisfied one multidegree at a time.

\subsection{Overview of algorithm}

The fundamental unit of computation is the routine \code{enumerate_tm} which
does this enlargement in a given multigraded component. More precisely, it takes
as input a representable family $I_t$ and a multidegree $(ijk)$ and produces a list of
representable families which together exhaust the set of ideals which (i)
are of the form $I_t + (p_1,\ldots,p_l)$ for some $t$ and $p_1,\ldots,p_l \in
S^i(A^*)\ot S^j(B^*) \ot S^k(C^*)$ and (ii) satisfy condition (\ref{codim}) in
multidegree $(ijk)$. Moreover, the families output by \code{enumerate_tm} are
parameterized combinatorially by the set of leading monomials of polynomials in
this multidegree. In principle, this routine already is sufficient to perform
full ideal enumeration, as one could fix a list of all the multidegrees and for
each $k$ compute all representable families satisfying (\ref{codim}) in the
first $k$ multidegrees with generators only in these multidegrees. Eventually,
by the Hilbert basis theorem,
% is there more here, some theorem about hilbert functions?
condition (\ref{codim}) will be satisfied in all codimensions. This stopping
condition can be checked computationally by computing the Hilbert series of the
common leading term ideal of the family. 
% this is a little imprecise, maybe I really have to pass to irreducible
% parameter spaces and generic leading term ideals to fully justify the stopping
% condition
In practice, we make the enumeration
more tractable via the use of two other routines \code{enumerate_tm_bounded} and
\code{ideal_sum}.

The routine \code{enumerate_tm_bounded} takes as input a multidegree $(ijk)$ and
a representable family $I_t$, say where $t$ ranges over $X$. It computes the
family restricted to the subvariety $Y\subset X$ on which $\operatorname{codim}
(I_t)_{ijk} \le r$. In other words, $t$ is restricted to the subvariety on which
a necessary condition holds for $I_t$ in the family to be contained in one
satisfying (\ref{codim}). We speak of \code{enumerate_tm_bounded} as applying
the $(ijk)$ test.
The routine \code{ideal_sum} takes a pair of families $I_t$ and $J_s$, say where
$t\in X$ and $s\in Y$, and computes the family $I_t + J_s$, where $(t,s) \in X
\times Y$.

A typical approach for border rank lower bounds is to enumerate representable
families satisfying (\ref{codim}) in multidegree (110) and to apply the (210)
and (120) tests to each family. Call the results of this the (110)
\emph{candidates}. For showing $2\times 2$ matrix multiplication has border rank
at least 7, this is already enough, as there are no candidates already at this
stage when $r=6$. If we need to continue, we similarly compute the (101) and
(011) candidates. For each triple of candidates in these multidegrees, we add
them together and apply the (111) test. This approach shows that $3\times 3$
matrix multiplication and the $3\times 3$ determinant polynomial each have 
border rank at least 17.

\subsection{Input description}\label{sec:input}
Choose bases of $A$, $B$ and $C$ consisting of weight vectors under the
torus of $B$. The data of the action of $B$ then consists of the weights of
these distinguished basis vectors along with a faithful representation $\mathfrak{n} \to
\mathfrak{gl}(A) \oplus \mathfrak{gl}(B)\oplus \mathfrak{gl}(C)\subset
\mathfrak{gl}(A\oplus B\oplus C)$ of the nilradical $\mathfrak{n}$ of the lie
algebra of $B$. Concretely, $\mathfrak{n}$ and its representation is described
by distinguishing lie algebra generators $x_1,\ldots,x_k \in \mathfrak{n}$ and
providing the matrices of the corresponding actions on the distinguished basis
of $A\oplus B\oplus C$. This data describing the action of $B$ along with $T$
expressed in the distinguished bases is the input to the routines in this file.
% TODO make code correspond to this paragraph. Namely, xs act on primal not on dual

\subsection{Representable families of ideals}\label{sec:representation}

\newcommand{\ba}{\textbf{a}} 
\newcommand{\bb}{\textbf{b}}
\newcommand{\bc}{\textbf{c}} 

We wish to describe a computationally efficient
representation of Borel fixed multigraded ideals which may depend on parameters.
At core, ideals will be
represented by Gr\"obner bases with homogeneous generators taken from
$T^\perp$. Thus, condition (\ref{inTperp}) will be satisfied by construction. 
The issues which must be addressed then are how to efficiently require ideals
are Borel fixed and how to represent parameter
dependence. 


\subsubsection{Borel fixed condition} 
\label{sec:borelfixed}
Immediately we can ensure ideals considered are torus invariant by
insisting generators are weight vectors. In other words the ideals under
consideration are homogeneous with respect to a finer grading determined by the
torus weights and can be taken to be generated homogeneously.

Now to enforce condition (\ref{Bfixed}), it is necessary and sufficient to
insist that the ideal is closed under the Lie algebra action of the generators
$x_1,\ldots,x_k \in \mathfrak{n}$. One could address this issue in a simple way:
Whenever a new ideal is constructed, recursively consider all pairs of current
ideal generators with Lie algebra generators $x_1,\ldots,x_k$ and adjoin as a
generator the corresponding raising. However, this approach is both slow and
error prone if one tries to make it efficient. We would like to arrange that
Gr\"obner basis computations in the underlying CAS automatically adjoin new
raisings of generators which may be needed. This would remove this burden from
us and also affords us the benefit of the significant engineering effort already
expended to do this kind of task quickly.

Conveniently, this is possible after making some observations. Let $\ba = \dim
A$, $\bb = \dim B$, $\bc = \dim C$, and write $a_1,\ldots, a_\ba \in A^*$,
$b_1,\ldots, b_\bb \in B^*$, $c_1,\ldots, c_\bc \in C^*$ for dual bases to the
distinguished bases of \S\ref{sec:input}. For notational convenience, also let
$n = \ba+\bb+\bc$ and for $1\le i\le n$, let $y_i$ range over the $n$ elements $
a_1,\ldots,a_\ba,b_1,\ldots,b_\bb,c_1,\ldots,c_\bc$. Let $x\in \mathfrak{n}$,
and suppose $x . y_j = \sum_{i=1}^\ba x_{ij} y_i$. The matrix $[x_{ij}]$ is
block diagonal with three diagonal blocks corresponding to the variables of
$A^*$, $B^*$, and $C^*$. Define the associated element $\bar x = \sum_{i,j=1}^n
x_{ij} y_i \partial_{y_j} $ in the Weyl algebra
$W=\CC[y_1,\ldots,y_n,\partial_{y_1},\ldots,\partial_{y_n}]$.

\begin{observation}\label{lieaction}
  Suppose $x\in \mathfrak{n}$ and $p \in \CC[y_1,\ldots,y_n] \subset W$.
  Then, in $W$, $\bar x p - p \bar x = x . p$.
\end{observation}
\begin{observation}\label{envelopingalgebra}
  If $x_1,\ldots,x_k$ generate $\mathfrak{n}$, then the subring of $W$ generated
  by $\bar x_1,\ldots, \bar x_k$ is the universal enveloping algebra
  $U(\mathfrak{n})$ 
\end{observation}

Let $S$ denote the subalgebra of $W$ generated by the $y_i$ and $\bar x_i$,
where $x_i$ generate $\mathfrak{n}$, and let $R = \CC[y_1,\ldots,y_n] \subset S$
denote the commutative polynomial ring. Then from Observation \ref{lieaction} it follows

\begin{observation}
  If $J \subset S$ is a two sided ideal, then $J\cap R$ is closed under the
  action of $\mathfrak{n}$. Moreover, if $I\subset R$ is an ideal then $SIS \cap R$ is the
  smallest $\mathfrak{n}$-fixed ideal containing $I$.
\end{observation}

In particular, to compute with only with Borel fixed multigraded ideals, it
suffices to work with two sided ideals of $S$ homogeneously generated by
elements of $R$. How do we work with such objects computationally? If now we
distinguish a vector space basis $x_1,\ldots, x_k$ of $\mathfrak{n}$, then $S$
is the noncommutative polynomial ring with indeterminants
$y_1,\ldots,y_n,x_1\ldots,x_k$ subject to the commutator relations $x_iy_j =
y_jx_i + x_i. y_j$, $x_ix_j = x_jx_i + [x_i,x_j]$, and the $y_i$ commute with
each other. This data makes $S$ into what is called a $G$-algebra
\cite{levandovsky}, and,
in particular, efficient implementations exist for computing Gr\"obner bases of
two sided ideals of $S$.

% There is an element in the Weyl algebra in the indeterminants
% $ a_1,\ldots,a_\ba,b_1,\ldots,b_\bb,c_1,\ldots,c_\bc$ corresponding to any
% element of $\mathfrak{n}$. Namely,

% $\mathbb{C}[a_1,\ldots,a_\ba,b_1,\ldots,b_\bb,c_1,\ldots,c_\bc,
% \partial_{a_1}, \ldots, \partial_{a_\ba},
% \partial_{b_1}, \ldots, \partial_{b_\bb},
% \partial_{c_1}, \ldots, \partial_{c_\bc} ]$ 

\subsubsection{Representation of parameterized families}
  \label{sec:parameterized}

Generators of a parameterized family $I_t \subset \CC[y_1,\ldots,y_n]$ are
essentially represnted as having coefficients in some other polynomial ring, and
we imagine the family to be parameterized by the different choices of values for
the new variables. In other words, a parameterized family $I_t$ is represented
by an ideal $J \subset \CC[y_1,\ldots,y_n,t_1,\ldots,t_k]$. We take the new
indeterminants $t_1,\ldots,t_k$ to have degree zero, so that a homogeneous
family $I_t$ is represented by a homogeneous ideal $J$.

For an family $I_t$ represented by $J$, if $t$ takes a value outside the variety
$X$ cut out by $J\cap \CC[t_1,\ldots,t_k]$, then by the Nullstellensatz $I_t =
(1)$. We therefore adopt the convention that the family $I_t$ corresponding to
$J$ takes parameter values $t$ only inside $X$.
In view of this convention, we say that a family $J'$ is a \emph{restriction} of
$J$ if, for $I' = J'\cap \CC[t_1,\ldots,t_k]$, we have $J' =
J+\CC[y_1,\ldots,y_n,t_1,\ldots,t_k]I'$. Hence, $J$ and $J'$ represent the same
family, but $J'$ is restricted to $X' = V(I') \subset X$.
For $t\in X$,
write $J(t)$ for the ideal in $\CC[y_1,\ldots,y_n]$ resulting from substituting
$t$ into each element of $J$, that is, $J(t)$ is the ideal $I_t$ of the
corresponding family. 

% TODO note semicontinuity of codim in t?

Fix a monomial ordering on $\CC[y_1,\ldots,y_n]$.
Condition (\ref{codim})
requires us to reason about Hilbert series of homogeneous ideals, and we will do
so by reducing to Hilbert series of the leading term ideal with respect to this
monomial ordering. 
We wish to choose a monomial
order on $\CC[y_1,\ldots,y_n,t_1,\ldots,t_k]$ so that the leading term ideal of
$J$ gives us information about the leading term ideals of each $J(t)$ in the
corresponding family.
It will be sufficient for our purposes to require that the monomial order is a block order
with the variable block $y_1,\ldots,y_n$ greater than the variable block
$t_1,\ldots,t_k$. In other words, we require that $y^{\alpha_1} t^{\beta_1} >
y^{\alpha_2} t^{\beta_2}$ if and only if $y^{\alpha_1} > y^{\alpha_2} $ or
$y^{\alpha_1} = y^{\alpha_2}$ and $t^{\beta_1} > t^{\beta_2}$ for some
fixed monomial order on just the $t$ monomials. 

This monomial order has the effect of grouping together monomials
with the same exponent in the variables $y_1,\ldots,y_n$, so it is possible to 
read off the coefficient of a leading monomial in the variables $y_1,\ldots,y_n$
as a polynomial in $\CC[t_1,\ldots,t_k]$. More precisely, for $p\in
\CC[y_1,\ldots,y_n,t_1,\ldots,t_k]$, if $\LM(p) = y^\alpha t^\beta$, write
$\Lm(p) = y^\alpha$, and denote by $\ltcoeff(p) \in
\CC[t_1,\ldots,t_k]$ the coefficient of $y^\alpha$ in $p$.

We need to understand how $\LT(J(t))$ can depend on $t$. Say that a
monomial $y^\alpha$ is \emph{trivially missing} if $y^\alpha t^\beta \in \LT(J)$
only if $t^\beta \in \LT(J)$. 
If $J\cap \CC[t_1,\ldots,t_k]$ is prime (maximal), for
instance, then $y^\alpha$ is trivially missing if and only if $y^\alpha \notin
\LT(J(t))$ for generic $t$ (for the unique $t$). 

% two concerns here (1) how to restrict t to that subvariety so that y^alpha not
% in LT(J(t)) (2) what happens under restriction


% need restriction to not  break test already done. actually not a soundness issue though
% (completeness)
% can always just apply necessary tests
% prove all ideals occur at some point in this process also need

We need the following:
\begin{observation}\label{missing}
  Let $J'$ be a restriction of $J$. % other conditions?
  If $s$ of the largest $m$ monomials in multidegree $(ijk)$ are trivially
  missing from $J$, then at least $s$ of the largest $m$ monomials are trivially
  missing from $J'$.

  In particular, if $s$ of the largest $m$ monomials are trivially missing in
  $J$, then at least $s$ of the largest $m$ monomials are not contained in
  $\LM(J(t))$ for any $t$.
\end{observation}
\begin{proof}
  Let $l = m-s$, $I = J\cap \CC[t_1,\ldots,t_k]$, and $I' = J'\cap
  \CC[t_1,\ldots,t_k]$.
  Suppose $p_1,\ldots,p_l\in J'$ have $\Lm(p_i)$ pairwise disjoint among the
  largest $m$ monomials of multidegree $(ijk)$ and $\ltcoeff(p_i) \notin I'$.
  It will be sufficient to show
  there exist $l$ polynomials satisfying the same conditions with $J'$ and $I'$
  replaced by $J$ and $I$, respectively. For each $i\le l$, there is $q_i$
  in the extension of $I'$ so that $p_i' = p_i+q_i \in J$. 
  We may assume $\ltcoeff(p_i') \notin I$.

  If $\Lm(p_i')$ are not pairwise distinct, suppose $\Lm(p_i') = \Lm(p_j')$ with
  $\Lm(p_i) > \Lm(p_j)$.
  Replace $p_i'$ with $p_i' \ltcoeff(p_j') - p_j' \ltcoeff(p_i')$.
  TODO maybe assume $I'$ maximal and argue we never get zero during this
  procedure
  % procedure assuming I' is maximal and considering smallest (I')^m
  % coefficients in

  % There are no $\CC[t_1,\ldots,t_k]$ syzygies between the $p_i'$. Indeed,
  % suppose $c_1p_1' + \cdots c_l p_l' = 0$, $c_i \in \CC[t_1,\ldots,t_k]$.
  % The coefficient of $\Lm(p_1)$ in each of $c_ip_i'$, $i > 1$ is contained in
  % $I'$, hence also $c_1p_1' \in I'$. Since $I'$ is prime, we have $c_1\in I'$.
  % We may argue similarly by induction that all $c_i \in I'$.
\end{proof}

Let $M$ be a set of monomials of multidegree $(ijk)$.
If a family $J$ contains a polynomial $p$ with $\LT(p) = y^\alpha$ for each
$y^\alpha\in M$ and all monomials of multidegree $(ijk)$ not in $M$
are trivially missing, then in view of observation \ref{missing}, all
ideals $J(t)$ have precisely the monomials of $M$ in the leading
term ideal in multidegree $(ijk)$. 
This fact is the basis of how \code{enumerate_tm} 
will enumerate all ideals of required codimension in a prescribed multidegree.
Central to this is the ability to restrict a given family to one where a
certain monomial is trivially missing.
If $y^\alpha$ is not trivially missing from $J$, then there is
$p\in J$ so that $\Lm(p) = y^\alpha$ and $\ltcoeff(p) \notin J\cap
\CC[t_1,\ldots,t_k]$. 
In particular, in order that a restriction $J'$ of $J$ have $y^\alpha$ trivially
missing, it is necessary for $\ltcoeff(p) \in J'$. Hence, the minimal
restriction $J'$ of $J$ where $y^\alpha$ is trivially missing can be computed by
recursively adjoining such leading terms to $J$. This is carried out in \code{minimal_relations_killing_leading_term}.



% Given a parameterized family $J(t)$ and a monomial $y^\alpha$, let $J_\alpha^0 =
% (\ltcoeff(p) \mid \Lm(p) \le \alpha \text{ and $p$ appears in the reduced
% Gr\"obner basis of $J$}) \subset \CC[t_1,\ldots,t_k]$. For instance, $J_0^0 =
% \CC[y_1,\ldots,y_n] \cap J$ is the set of defining equations of the domain of
% $t$.
% It is clear that $t \in V(J_\alpha^0)$ is a necessary condition that $y^\alpha
% \notin \LT(J(t))$. It need not be sufficient, as under some parameter values,
% $\ltcoeff$'s of polynomials in the ideal may vanish, promoting other monomials
% dividing $y^\alpha$ to be leading. % this might be misleading

% Write $J_\alpha = J+J_\alpha^0$, which corresponds to restricting the family $J(t)$
% to $t\in V(J_\alpha^0)$. We can iterate this restriction process until
% stabilization to obtain $J_\alpha^* = J_{\alpha\cdots\alpha}$.

% However, we have the following, which will be sufficient for our purposes.

% \begin{observation}
%   Suppose $y^\alpha$ is a monomial satisfying the property that $y^\alpha
%   t^\beta \in \LT(J) \implies t^\beta \in \LT(J)$ for every $\beta$. Then
%   $y^\alpha \notin \LT(J(t))$ for generic $t$.
% \end{observation}


% \begin{observation}
%   Let $y^\alpha$ be a monomial. If, for each monomial $y^\beta$ 
%   of multidegree strictly smaller than that of $y^\alpha$, $\ltcoeff(y^\beta)$
%   is either $J_0$ or $(1)$, then 

  
% \end{observation}

% The hypotheses of the observation can be interpreted as saying $\LT(J)$ is
% independent of parameters inside multigraded components

% With this choice of monomial ordering, some facts will be useful.
% Let $\Sigma$ be an order ideal of multidegrees, that is a set 
% closed under taking smaller multidegrees, and let $R_\Sigma$ be the sum of the 
% multigraded components corresponding to $\Sigma$.
% We say a family $J(t)$ has
% \emph{leading term ideal independent of parameters on $\Sigma$} if $\LT(J(t))
% \cap R_\Sigma$ is independent of the choice of $t$. 

% \begin{observation}
%   $J(t)$ has leading term ideal independent of parameters on $\Sigma$
%   if and only if $\LT(J)$ has a set of generators $y^{\alpha_i} t^{\beta_i}$
%   where $\beta_i = 0$ when $\alpha_i \in \Sigma$.
%   % In this case, $\LT(J(t)) \cap R_\Sigma = \LT(J)\cap R_\Sigma$.
% \end{observation}

% \begin{observation}
%   Let $J\cap \CC[t_1,\ldots,t_k]$ be prime. Then $\LT(J(t))$ for a generic $t$
%   depends only on $\LT(J)$. 
% \end{observation}

% \begin{observation}
%   Let $y^\alpha$ be a monomial and let $\Sigma$ be the set of multidegrees
%   stictly smaller than the multidegree of $y^\alpha$.
%   If $J$ has leading term ideal independent of parameters on $\Sigma$, 

%   Let $p$ be an element of the reduced Gr\"obner basis of $J$, and $\Sigma$ the
%   set of multidegrees strictly smaller than the degree of $p$. 
% \end{observation}

% % For homogeneous ideals, this property of $I_t$ only depends on the elements of
% % $I_t$ contained in those components


% % Over the course of the ideal enumeration algorithm, it will be required to
% % compute the parameter values $t$ for which a given monomial does not occur in
% % the leading term ideal $I_t$ . 

% % \begin{observation}
% % If
% %   Suppose $X$ is irreducible. Then 
% %   $m t^\alpha \in \LT(J)$ for some $\alpha$ if and only if
% %   $m \in \LT(J(t))$ for generic $t\in X$.
% % \end{observation}
% % \begin{observation}
% %   Suppose $X$ is irreducible. Then 
% %   $m t^\alpha \in \operatorname{LT}(J)$ for some $\alpha$ if and only if
% %   $m \in \operatorname{LT}(J(t))$ for generic $t\in X$.
% % \end{observation}

% % TODO Are there natural conditions on J such that if $X$ is irreducible then
% % the family has constant hilbert function on all of X or just an open subset?

% This observation suggests an algorithm for the task.

% (see, e.g.,
% \code{minimal_relations_killing_leading_term})

\subsubsection{Putting it together}\label{sec:representationlast}

To combine the ideas of \S \ref{sec:borelfixed} and \S \ref{sec:parameterized}
we take as our working ring
$S=\CC[y_1,\ldots,y_n,x_1,\ldots,x_k,t_1,\ldots,t_{k'}]$, where $y_i$ denote the
variables of $A^*$, $B^*$, and $C^*$, $x_i$ are identified with a vector space
basis of $\mathfrak{n}$, and $t_i$ are extra parameters of degree and weight 0.
Here the $y_i$ and $x_i$ satisfy the noncommutative relations described in \S
\ref{sec:borelfixed}, and both sets commute with the $t_i$. The block condition
on the monomial order of \S\ref{sec:borelfixed} here requires that the variable
block $y_1,\ldots,y_n,x_1,\ldots,x_k$ is greater than the variable block
$t_1,\ldots,t_{k'}$.
A family of ideals $I_t$ is representable (see
\S\ref{sec:intro}) when it corresponds to a two-sided ideal $J\subset S$
homogeneously generated (including weight homogeneously) by elements of the
subring $\CC[y_1,\ldots,y_n,t_1,\ldots,t_{k'}]$. 
% TODO technically to agree with intro I should require $J \subset Tperp
The notions of restriction of
families, evaluation at $t$, and trivially missing monomials are carried over
unmodified from \S\ref{sec:parameterized}.

In multidegree $(ijk)$ pick an ordering $(m_i)_i$ of the monomials $m_i =
y^{\alpha_i}$ so that if $\wt(y^{\alpha_i}) \le \wt(y^{\alpha_j})$ in the
partial order, then $i \ge j$ (note the order reversal), and if
$\wt(y^{\alpha_i}) = \wt(y^{\alpha_j})$ and $y^{\alpha_i} < y^{\alpha_j}$,
then $i \le j$. Hence, the last monomial on this list is the biggest monomial of
the smallest weight.
Note, there is no harm in choosing the monomial order on $y_1,\ldots,y_n$ so
that it already satisfies these properties. However, this is not strictly
necessary and could impact the performance of the implementation, so we do not
make this requirement.

This order has the property that
adjoining a homogeneous polynomial $p$ to an
ideal $J$ doesn't change the sets of leading coefficients corresponding to
larger monomials than $\Lm(p)$ in $(m_i)_i$.
That is, for $p$ with $\Lm(p) = m_i$ and $j>i$, the sets
$\{\ltcoeff(q) : \Lm(q) = m_j\}$ corresponding to $J$ and $J+(p)$ are 
identical. 
For such an ordering, and under our new conditions on $J$ we have the following
analogue of observation \ref{missing}:

\begin{observation}\label{missing2}
  Let $J'$ be a restriction of a representable family $J$. % other conditions?
  If $s$ of the largest $m$ monomials of $(m_i)_i$ in multidegree $(ijk)$ are trivially
  missing from $J$, then at least $s$ are trivially missing from $J'$.

%   In particular, if $s$ of the largest $m$ monomials are trivially missing in
%   $J$, then at least $s$ of the largest $m$ monomials are not contained in
%   $\LM(J(t))$ for any $t$.
\end{observation}

The routine \code{get_vorder} chooses a weight basis of $(T^\perp)_{ijk}$ with
pairwise distinct leading monomials and orders them according to the
requirements for the $m_i$ (thus implicitly fixing such a list $(m_i)_i$). 
These orderings are used to organize the calculations in \code{enumerate_tm} and
\code{enumerate_tm_bounded}.

\subsection{Software used}

The code is written in Python 3, using the libraries provided by SageMath.
The implementation of $G$-algebras and their arithmetic and ideal computations
is done under the hood by Plural, a kernel extension of the Singular computer
algebra system.

\section{Program Listing}
"

<<complete=False>>=
from itertools import *
from operator import sub
from functools import cache

@

<<complete=False>>=

def lp_look_strict_by_deg(ba,deg_bound=2,Jxi=None):
    from operator import concat
    from sage.numerical.mip import MIPSolverException
    def dfs(Jxi,d):
        if d > deg_bound:
            yield Jxi
            return
        J,S,xi = Jxi
        core = ba.ideal_constant_generators(J)
        mod = ba.zero if core.is_zero() else JDat(ba,core) 
        tms = tuple(sorted(map(tuple,IntegerVectors(d, len(ba.dims)))))
        tmcs = tuple(sorted(map(tuple,IntegerVectors(d+1, len(ba.dims)))))
        for Jxiup in ba.get_ideals_for_tm(tms, tmcs, J, S, xi, 
                lp_interval_length_bound=max(3, 7-d), precise_lp=True, mod=mod):
            for res in dfs(Jxiup,d+1):
                yield res
    return dfs((ba.R0.ideal(side='twosided'),[ba.R0.one()], 0) if Jxi is None else Jxi,1)

def lp_look_strict(ba,deg_bound=2,tm_max_bound=oo,Jxi=None):
    from operator import concat
    from sage.numerical.mip import MIPSolverException
    @cache
    def get_tms(d,tmma):
        tms = sorted(map(tuple,IntegerVectors(d, len(ba.dims))))
        return tuple([tm for tm in tms if max(tm) == tmma])
    def dfs(Jxi,d,tmma):
        if d > deg_bound:
            yield Jxi
            return
        J,S,xi = Jxi
        core = ba.ideal_constant_generators(J)
        mod = ba.zero if core.is_zero() else JDat(ba,core) 
        tms = get_tms(d,tmma)
        tmcs = tuple((tm for tma in range(1,tmma+2) for tm in get_tms(d+1,tma)))
        for Jxiup in ba.get_ideals_for_tm(tms, tmcs, J, S, xi, 
                lp_interval_length_bound=max(3, 7-d), precise_lp=False, mod=mod):
            for res in dfs(Jxiup,d+1 if tmma == min(d,tm_max_bound) else d, 
                           1 if tmma == min(d,tm_max_bound) else tmma+1):
                yield res
    return dfs((ba.R0.ideal(side='twosided'),[ba.R0.one()], 0) if Jxi is None else Jxi,1,1)

def three_place_ideals_up_to(ba,d=1,Jxi=None):
    def dfs(Jxi,tms,mod=ba.zero): 
        if len(tms) == 0:
            yield Jxi
            return
        tm = tms[0]
        print (tm)
        J,S,xi = Jxi
        for Jxiup in ba.get_ideals_for_tm([tm], [tm[:s]+(tm[s]+1,)+tm[s+1:] for s in range(3)]
                      +([(2,2,2)] if sum(tm) == 4 and all(e <= 2 for e in tm) else []),
                  J,S,xi,search_above = sum(tm) == 2, lp_interval_length_bound = 6
                                          if sum(tm) <= 3 else 2, mod = mod):
            if len(tms) == 1 and sum(tm) < d:
                ltIlb, ltIub = ba.leading_term_ideal(*Jxiup[:2])
                R = ltIlb.ring()
                
                degseries = ltIlb.hilbert_series() == ba.target_hilbert_series_td() and \
                    ltIub.hilbert_series() == ba.target_hilbert_series_td()
                
                for dnext in range(sum(tm)+1,d+1):
                    print (dnext)
                    tmsnext = []
                    if dnext > 3 and degseries:
                    # if dnext > sum(tm)+1 and dnext > 3 and degseries:
                        break
                    for tmc in IntegerVectors(dnext,3):
                        tmc = tuple(tmc)
                        codim_upper_bound = len([m for m in mod.get_vs(tmc) if liftR(m.lm(),R) not in ltIlb]) + ba.tm_codim_tperp(tmc)
                        codim_noadd_lower_bound = len([m for m in mod.get_vs(tmc) if liftR(m.lm(),R) not in ltIub]) + ba.tm_codim_tperp(tmc)
                        if codim_upper_bound != ba.r or codim_noadd_lower_bound != ba.r:
                            tmsnext.append((tmc,codim_upper_bound))
                    if len(tmsnext) > 0:
                        break
                tmsnext.sort(key = lambda p: (p[1], max(p[0]), p[0]))
                print(tmsnext)
                tmsnext = [tmc for tmc,_ in tmsnext]
                modnext = JDat(ba,ba.ideal_constant_generators(Jxiup[0]))
            else:
                tmsnext = tms[1:]
                modnext = mod
                
            for Jxifull in dfs(Jxiup,tmsnext,modnext):
                    yield Jxifull

    if Jxi is not None:
        mod = JDat(ba,ba.ideal_constant_generators(Jxi[0]))
        for Jxifull in dfs(Jxi, [(1,0,0),(0,1,0),(0,0,1)], mod):
            yield Jxifull
        return

    cd2 = [ ba.get_ideals_for_tm([(1,1,0)],tmcs=[(2,1,0),(1,2,0)]) ]
    cd2.append(ba.get_ideals_for_tm([(1,0,1)],tmcs=[(2,0,1),(1,0,2)]))
    cd2.append(ba.get_ideals_for_tm([(0,1,1)],tmcs=[(0,2,1),(0,1,2)]))

    for Jxis in lproduct(*map(enumerate,cd2)):
        print ('triple',[i for i,Iyi in Jxis])
        Jxi = ba.sum_ideals_S([Iyi for i,Iyi in Jxis])
        try:
            next(ba.get_ideals_for_tm([],[(1,1,1)],*Jxi))
            for Jxifull in dfs(Jxi,[(2,0,0),(0,2,0),(0,0,2)] if d >= 2 else []):
                yield Jxifull
        except StopIteration:
            pass
        
def lp_only_look_nonstrict(ba):
    from operator import concat
    from sage.numerical.mip import MIPSolverException
    tms = [(1,1,0),(1,0,1),(0,1,1)]
    tmcss = [((2,1,0),(1,2,0)), 
              ((2,0,1),(1,0,2)),
              ((0,2,1),(0,1,2))]
    ips = [ IntegerProgram(ba.zero.get_linear_program((tm,),tmcs,interval_length_bound=8),
                           [p.lm() for p in ba.zero.get_vs(tm)])
            for lpi,(tm,tmcs) in enumerate(zip(tms,tmcss)) ]
    def prunef(psol):
        psols = [[] for _ in tms]
        for x,v in psol.items():
            psols[tms.index(ba.polynomial_tm(x))].append((x,v))
        for psol, ip in zip(psols,ips):
            if not ip.can_extend_psol(psol):
                return False
        return True
    lp = ba.zero.get_linear_program(((1,1,0),(1,0,1),(0,1,1)),((1,1,1),),
                                    interval_length_bound=4)
    xs = [x.lm() for tm in tms for vs in ba.zero.get_vss(tm).values() for x in vs]
    
    @cache
    def getJxis(lpi,sol):
        tm = tms[lpi]
        tmcs = tmcss[lpi]
        return list(ba.get_ideals_for_tm((tm,),tmcs,preset_choices=dict(sol)))
    
    for fsol in lp_integer_points(lp,xs,fullsol=False,prunef=prunef):
        sols = [{} for _ in tms]
        for x,v in fsol.items():
            sols[tms.index(ba.polynomial_tm(x))][x] = v
        Jxiss = []
        for Jxis in product(*[getJxis(lpi,tuple(sorted(sol.items()))) 
                              for lpi, sol in enumerate(sols)]):
            Jxi = ba.sum_ideals_S(Jxis)
            try:
                next(ba.get_ideals_for_tm([],[(1,1,1)],*Jxi))
                yield Jxi
            except StopIteration:
                pass

def lp_only_look_strict(ba):
    from operator import concat
    from sage.numerical.mip import MIPSolverException
    tms = [(1,1,0),(1,0,1),(0,1,1)]
    tmcss = [((2,1,0),(1,2,0)), 
              ((2,0,1),(1,0,2)),
              ((0,2,1),(0,1,2))]
    lp = deepcopy(ba.zero.get_linear_program(
        ((1,1,0),(1,0,1),(0,1,1)),((1,1,1),),interval_length_bound=4))
    sols_ss = []
    for tm,tmcs in zip(tms,tmcss):
        lpcur = ba.zero.get_linear_program((tm,),tmcs,interval_length_bound=6)
        xs = [x.lm() for vs in ba.zero.get_vss(tm).values() for x in vs]
        ss = SetSystem()
        sols_ss.append(ss)
        R = PolynomialRing(QQ,'x',len(xs))
        bad_condition = R.ideal(1)
        for si,sol in enumerate(lp_integer_points(lpcur,xs,fullsol=False)):
            print(tm,si)
            sol = [m for m,v in sol.items() if v==0]
            ss.add_set(sol)
            bad_condition = bad_condition.intersection(R.ideal([
                R.gen(xs.index(m)) for m in sol]))
        print(bad_condition)
        for bad in bad_condition.gens():
            assert len(bad.exponents()) == 1
            assert all(k==1 for i,k in bad.exponents()[0].sparse_iter())
            lp.add_constraint(lp.sum(lp[xs[i].lm()] for i,_ in bad.exponents()[0].sparse_iter()) <= bad.degree()-1)
    @cache
    def getJxis(lpi,sol):
        tm = tms[lpi]
        tmcs = tmcss[lpi]
        preset_choices = {v.lm() : True for v in ba.zero.get_vs(tm)}
        preset_choices.update({m.lm() : False for m in sol})
        res = list(ba.get_ideals_for_tm((tm,),tmcs,preset_choices=preset_choices))
        if len(res) == 0:
            sols_ss[lpi].remove_set(sol)
        return res
    def prunef(psol):
        Slohis = [([],set([v.lm() for vs in ba.zero.get_vss(tm).values() 
                       for v in vs])) for tm in tms]
        for m,v in psol.items():
            lpj = tms.index(ba.polynomial_tm(m))
            if v == 0:
                Slohis[lpj][0].append(m)
            else:
                Slohis[lpj][1].remove(m)
        return all( any(True for _,si in ss.iter_sets(Slo,Shi) )
                   for ss,(Slo,Shi)
                   in zip(sols_ss,Slohis) )
    xs = [x.lm() for tm in tms for vs in ba.zero.get_vss(tm).values() for x in vs]
    for fsol in lp_integer_points(lp,xs,fullsol=False,prunef=prunef):
        sols = [[] for _ in tms]
        for x,k in fsol.items():
            if k == 0:
                sols[tms.index(ba.polynomial_tm(x))].append(x)
        sols = [tuple(sorted(sol)) for sol in sols]
        Jxiss = [getJxis(lpj,sol) for lpj,sol in enumerate(sols)]
        for Jxis in product(*Jxiss):
            Jxi = ba.sum_ideals_S(Jxis)
            try:
                next(ba.get_ideals_for_tm([],[(1,1,1)],*Jxi))
                yield Jxi
            except StopIteration:
                pass

def lp_only_look(ba):
    from operator import concat
    from sage.numerical.mip import MIPSolverException
    tms = [(1,1,0),(1,0,1),(0,1,1)]
    tmcss = [((2,1,0),(1,2,0)), 
              ((2,0,1),(1,0,2)),
              ((0,2,1),(0,1,2))]
    lps = []
    xss = []
    for tm,tmcs in zip(tms,tmcss):
        lp = ba.zero.get_linear_program((tm,),tmcs,interval_length_bound=7)
        # lp = ba.zero.get_linear_program(tuple(tms),tmcs+((1,1,1),))
        # lp = ba.zero.get_linear_program((tm,),tmcs+((1,1,1),))
        xs = [x.lm() for vs in ba.zero.get_vss(tm).values() for x in vs]
        lps.append(lp_integer_points(lp,xs,fullsol=False))
        xss.append(xs)
    lp = ba.zero.get_linear_program(((1,1,0),(1,0,1),(0,1,1)),((1,1,1),),interval_length_bound=4)
    xs = [x.lm() for tm in tms for vs in ba.zero.get_vss(tm).values() for x in vs]
    yss = [[] for _ in range(3)]
    sols = [[] for _ in range(3)]
    sols_ss = [SetSystem() for _ in range(3)]
    remove = set()
    @cache
    def getJxis(lpi,si):
        tm = tms[lpi]
        tmcs = tmcss[lpi]
        res = list(ba.get_ideals_for_tm((tm,),tmcs,preset_choices=sols[lpi][si]))
        if len(res) == 0:
            remove.add((lpi,si))
            print("REMOVE",remove)
        return res
    si = 0 
    while True:
        fails = 0
        for lpi,(lpx,xsx) in enumerate(zip(lps,xss)):
            try:
                sol = next(lpx)
                print(lpi,si)
                yss[lpi].append((lpi,si))
                sols[lpi].append(sol)
                def getSlohis(psol):
                    Slohis = [([],set([v.lm() for vs in ba.zero.get_vss(tm).values() 
                                   for v in vs])) for tm in tms]
                    for m,v in psol.items():
                        lpj = tms.index(ba.polynomial_tm(m))
                        if v == 0:
                            Slohis[lpj][0].append(m)
                        else:
                            Slohis[lpj][1].remove(m)
                    return Slohis
                def prunef(psol):
                    return all( any(True for _,si in ss.iter_sets(Slo,Shi) 
                                     if (lpj,si) not in remove)
                               for lpj,(ss,(Slo,Shi))
                               in enumerate(zip(sols_ss,getSlohis(psol))) if lpi != lpj)
                for m,v in sol.items():
                    if v == 0:
                        lp.set_max(lp[m],0)
                    else:
                        lp.set_min(lp[m],1)
                try:
                    lp.solve()
                    sols_ss[lpi].add_set([m for m,v in sol.items() if v==0],si)
                    for fsol in lp_integer_points(lp,xs,fullsol=False,prunef=prunef):
                        Slohis = getSlohis(fsol)
                        # assert all(sorted(Slo) == sorted(Shi) for Slo,Shi in Slohis)
                        sis = [next(ss.iter_sets(Slo,Shi))[1]
                               for ss,(Slo,Shi) in zip(sols_ss,Slohis)]
                        Jxiss = [getJxis(lpj,si) for lpj,si in enumerate(sis)]
                        assert len(Jxiss) == 3
                        for Jxis in product(*Jxiss):
                            Jxi = ba.sum_ideals_S(Jxis)
                            try:
                                next(ba.get_ideals_for_tm([],[(1,1,1)],*Jxi))
                                yield Jxi
                            except StopIteration:
                                pass
                except MIPSolverException:
                    print("SKIPPING")
                    pass
                for m in sol.keys():
                    lp.set_max(lp[m],1)
                    lp.set_min(lp[m],0)
            except StopIteration:
                fails += 1
        if fails == 3:
            break
        si += 1
        
@

The Python class \code{BorderApolarity} provides a convenient interface to the
routines of this file which depend on the input data (see \S\ref{sec:input}).

<<complete=False>>=

class BorderApolarity:
    def __init__(self, Tinfo, r, basis=True):
        self.T, self.dims, self.vwts, self.xs = Tinfo
        self.xs = [-x.T for x in self.xs]
        self.vwts = -self.vwts
        if basis:
            self.Bxs = self.xs
        else:
            self.Bxs = basis_of_lie_algebra(self.xs)
        self.r = r
        self.F = self.T.base_ring()
        self.nvars = sum(self.dims)
        self.nvarsx = self.nvars + len(self.Bxs)
        self.R0 = self.get_ring()
        self.vixs = [0]
        for s in self.dims:
            self.vixs.append(self.vixs[-1] + s)
        self.Tdeg = [0 for i in range(len(self.dims))]
        e = self.T.exponents()[0]
        for i,k in e.sparse_iter():
            self.Tdeg[ next(j for j in range(len(self.dims)) 
                            if i < sum(self.dims[:j+1])) ] += k
        self.Tdeg = tuple(self.Tdeg)

        self.zero = JDat(self,self.R0.ideal())
        self.get_ring_exact = cache(self.get_ring_exact)
        self.tm_dim = cache(self.tm_dim)
        self.tm_dim_tperp = cache(self.tm_dim_tperp)
        self.tm_codim_tperp = cache(self.tm_codim_tperp)
        self.tm_target_out = cache(self.tm_target_out)
        self.target_hilbert_series = cache(self.target_hilbert_series)
        self.target_hilbert_series_td = cache(self.target_hilbert_series_td)
        self.Tperp = cache(self.Tperp)
    @

The function \code{get_ring} implements the logic described in
\S\ref{sec:representation}, namely, it forms the $G$-algebra generated by the
variables $a_1,\ldots,a_\ba$, $b_1,\ldots,b_\bb$,
$c_1,\ldots,c_\bc$, $x_1,\ldots,x_k$, $t_1,\ldots,t_{\text{params}}$ subject to the
appropriate commutator relations and with 
degree reverse lexicographic monomial order on the variables 
$a_1,\ldots,a_\ba$, $b_1,\ldots,b_\bb$,
$c_1,\ldots,c_\bc$, $x_1,\ldots,x_k$ and $t_1,\ldots,t_{\text{params}}$ seperately,
with products of monomials between these groups ordered lexicographically. 

The computation of this algebra is expensive, so 
for performance reasons, calling this function actually chooses the
next power of two number of parameters larger than the requested number and
caches the results for the future. 

<<complete=False>>=
    def get_ring(self,params=0):
        if params > 0:
            params = next(1<<k for k in range(32) if 1<<k >= params )
        return self.get_ring_exact(params)

    def get_ring_exact(self, params=0):
        print("get_ring_exact",params)
        if params > 0:
            params = next(1<<k for k in range(32) if 1<<k >= params )
        # should give cached rings of power of two size
        lets = 'abcdefghijklmnopqrsuvwyz'
        B = matrix([x.list() for x in self.Bxs]).T
        A = FreeAlgebra(self.F,['%s%d' % (lets[i],j) for i in range(len(self.dims))
            for j in range(self.dims[i])] + 
                ['x%d'%i for i in range(len(self.Bxs))] + ['t%d' % j for j in range(params)])
        yvs = A.gens()[:sum(self.dims)]
        xvs = A.gens()[sum(self.dims):sum(self.dims)+len(self.Bxs)]
        rels = {}
        for (xv1,x1),(xv2,x2) in combinations(zip(xvs,self.Bxs),2):
            x3 = x2*x1 - x1*x2
            if not x3.is_zero():
                xv3 = sum(e*yv for e,yv in zip(B.solve_right(vector(x3.list())).list(),xvs))
                rels[xv2*xv1] = xv1*xv2 + xv3
        for (xv,x),(yi,yv) in product(zip(xvs,self.Bxs),enumerate(yvs)):
            xy = x*vector(self.F,self.Bxs[0].nrows(),{yi:1})
            if not xy.is_zero():
                xyv = sum(e*yv for e,yv in zip(xy.list(),yvs))
                rels[xv*yv] = yv*xv + xyv
        order = TermOrder('degrevlex',sum(self.dims)+len(self.Bxs))
        if params > 0:
            order = order + TermOrder('degrevlex',params)
        return A.g_algebra(rels,order=order)
@

\code{get_vorder} takes a multidegree $ijk$ as argument and 
computes a basis of 
$(T^\perp)_{ijk}$ as described at the end of \S\ref{sec:representationlast}.

<<complete=False>>=
    
    def x_weight(self, x):
        i,j = x.nonzero_positions()[0]
        wt = tuple(self.vwts.column(i)-self.vwts.column(j))
        for i,j in x.nonzero_positions()[1:]:
            assert wt == tuple(self.vwts.column(i)-self.vwts.column(j))
        return wt

    def weight_le(self, wta, wtb):
        dwt = vector(self.F,map(sub,wtb,wta))
        u = dwt * self.Xi
        return all(e >= 0 for e in u) and u*self.X == dwt

    def lm(self, p, toR0 = False):
        R = self.R0 if toR0 else p.parent()
        if p.is_zero():
            return R.zero()
        return R.prod(R.gen(x)**k for x,k in 
                p.lm().exponents()[0][:self.nvars].sparse_iter())
@


<<complete=False>>=
                   
    # adjoining these to the ideal will not necessarily zero out the
    # corresponding coefficient unless larger monomial
    # coefficients of the same weight and degree are not
    # affected
    # J assumed to be twostd()
    def leading_term_coefficients(self,J,p):
        ts = []
        for q in J.gens():
            if not self.lm(q).is_constant() and \
                    self.R0.monomial_divides(self.lm(q,True),self.lm(p,True)):
                ts.append(q.coefficient(self.lm(q)))
        return ts

    def get_param_ideal(self,J,xi=None):
        if xi is None:
            xi = J.ring().ngens() - self.nvarsx
        if J.ring().ngens() == self.nvarsx:
            assert xi == 0
            return None,None
        R = PolynomialRing(ba.F,J.ring().gens()[self.nvarsx:self.nvarsx+xi],implementation='singular')
        I = R.ideal([R(p) for p in J.gens() if self.lm(p).is_constant()])
        def coerce(t):
            return R(t)
        return I,coerce
    
    def family_dimension(self,J,xi):
        I,_ = self.get_param_ideal(J,xi)
        if I is None:
            return 0
        return I.dimension()

    def lift_param(self,R,t):
        tg = t.parent().ngens()
        Rg = R.ngens()-self.nvarsx
        assert Rg >= tg or all(next(tt.exponents()[0].sparse_iter())[0] < Rg for tt in t.variables())
        return R(t(R.gens()[self.nvarsx:self.nvarsx+tg]+(R.zero(),)*(tg-Rg)))
    
    def get_ideals_for_tm(self, tms=[], tmcs=[], J=None, S=None, xi=0,
                          search_above=True, lp_interval_length_bound=6, precise_lp=False,
                          preset_choices = {}, mod=None, sound=False):
        if mod is None:
            mod = self.zero
        tms = tuple(tms)
        tmcs = tuple(tmcs)
        
        tmmax = reduce(lambda a,b: tuple(map(max,a,b)),tms+tmcs)
        
        gbbound = max([sum(tm) for tm in tms+tmcs])
        precise_gb_bound = False
        sound_gb = sound or precise_gb_bound
        cur_search_above = search_above and sound_gb
        # searching above only sound if sound_gb, otherwise, might fail to
        # search possibility of including a leading monomial which is possibly
        # present but not seen in our generating set

        def gb(J):
            if not sound_gb or J.ring() is self.R0:
                return self.bounded_gb(J,gbbound)
            else:
                return J.twostd()
        
        if J is None:
            J = self.get_ring().ideal(side='twosided')
            S = [J.ring().one()]
            
        Jorig = J
        xiorig = xi
        if precise_gb_bound:
            J += [p for i,e in enumerate(tmmax) for p in
                    self.all_monomials((0,)*i+(e+1,)+(0,)*(len(tmmax)-i-1),J.ring())]
            
        from copy import deepcopy
        if sound:
            vkstart,vkbelowcnt,wkabove = mod.get_poset_info_for_search(
                    (), tms+tmcs if cur_search_above else tms)
        else:
            vkstart,vkbelowcnt,wkabove = mod.get_poset_info_for_search(
                    tms, tmcs if cur_search_above else ())
        vkstart = deepcopy(vkstart)
        vkbelowcnt = deepcopy(vkbelowcnt)
        tmsvks = set((tm,wt,i) for tm in tms for wt,vs in mod.get_vss(tm).items() 
                     for i in range(len(vs)))
        
        def guaranteed_lms(J,S):
            lts = SetSystem()
            if J.ring().ngens() == self.nvarsx:
                for p in J.gens():
                    if not p.is_zero():
                        lts.add_set([xi for xi,e in p.lm().
                                 exponents()[0].sparse_iter() for _ in range(e)],None)
                lms = [p.lm() for p in J.gens() if not p.is_zero()]
            else:
                I,coerce = self.get_param_ideal(J)
                SI = [coerce(s) for s in S]
                lms = []
                for p in J.gens():
                    if p.is_zero():
                        continue
                    It = I + coerce(p.coefficient(self.lm(p)))
                    if any(s in It for s in SI):
                        lts.add_set([xi for xi,e in self.lm(p,True).
                                 exponents()[0].sparse_iter() for _ in range(e)],None)
            def in_lti(m):
                return any(True for _ in lts.iter_sets([],[xi 
                       for xi,e in m.exponents()[0].sparse_iter() for _ in range(e)]))
            return in_lti
                        
        bp = BinaryProgram(mod.get_linear_program(tms,tmcs,
             interval_length_bound=lp_interval_length_bound,precise=precise_lp),
                           xs = [p.lm() for tm in tms for p in mod.get_vs(tm)],
                            use_infeas=False)
        lpset = set()
        unset = set()
        params_needed = {vks : vks[2] for vks in tmsvks}
        num_in = {tm:0 for tm in tms+tmcs}
        num_out = {tm:0 for tm in tms+tmcs}
        in_lti = guaranteed_lms(J,S)
        for tm in tms+tmcs:
            for wt,vs in mod.get_vss(tm).items():
                for i,p in enumerate(vs):
                    if in_lti(p.lm()):
                        bp.set_min(p.lm(),1)
                        lpset.add((p.lm(),1))
                        num_in[tm] += 1
                        if tm in tms:
                            for j in range(i+1,len(vs)):
                                params_needed[(tm,wt,j)] -= 1
                    else:
                        unset.add((tm,wt,i))
                        
        if bp.extend() is None:
            print ("get_ideals_for_tm: no solution to start")
            return
        
        J = gb(J)
        in_lti = guaranteed_lms(J,S)
        sav = []
        for vk in unset:
            tm,wt,i = vk
            vs = mod.get_vss(tm)[wt]
            p = vs[i].lm()
            if in_lti(p):
                sav.append(vk)
                bp.set_min(p.lm(),1)
                lpset.add((p.lm(),1))
                num_in[tm] += 1
                if tm in tms:
                    for j in range(i+1,len(vs)):
                        params_needed[(tm,wt,j)] -= 1
        for vk in sav:
            unset.remove(vk)
                        
        lpsols = SetSystem()
        nlpsols = 0
        tm_xs = [p.lm() for tm in tms for p in mod.get_vs(tm)]

        if bp.extend() is None:
            print ("get_ideals_for_tm: no solution to start needed gb")
            return
        
        processed = [0]
        height = [0]
        vstr_width = [0]
        def dfs(J,S,xi):
            # global s1
            # found = all((lp.get_min(lp[w.lm()]) == 1) ==
            #             (s1.reduce(w.lm()) != w.lm()) for w in
            #             mod.get_vs((1,1)) if lp.get_min(lp[w.lm()]) ==
            #             lp.get_max(lp[w.lm()]))
            if all(num_in[tm] == len(mod.get_vs(tm)) - self.tm_target_out(tm) and \
                           num_out[tm] == self.tm_target_out(tm) for tm in tms) and \
                       (not cur_search_above or all(num_out[tm] >= self.tm_target_out(tm) for tm in tmcs)):
                yield (J,S,xi)
                return

            def vkremove(vk,added=None):
                # print '%sremove1'%(' '*processed[0]),vk
                vkstart.remove(vk)
                for wk in wkabove[vk]:
                    vkbelowcnt[wk] -= 1
                    if vkbelowcnt[wk] == 0:
                        # print '%sadd2'%(' '*processed[0]),wk
                        vkstart.add(wk)
                        if added is not None:
                            added.append(wk)
                processed[0] += 1
                return added
            def vkadd(vk,removed=None):
                # print '%sadd3'%(' '*processed[0]),vk
                for wk in wkabove[vk]:
                    if vkbelowcnt[wk] == 0:
                        # print '%sremove2'%(' '*processed[0]),wk
                        vkstart.remove(wk)
                        if removed is not None:
                            removed.append(wk)
                    vkbelowcnt[wk] += 1
                # print '%sadd1'%(' '*processed[0]),vk
                vkstart.add(vk)
                processed[0] -= 1
                return removed
                    
            set_stack = []
            set_to_process = [vk for vk in vkstart if vk not in unset]
            while len(set_to_process) > 0:
                vk = set_to_process.pop()
                set_stack.append(vk)
                added = vkremove(vk,[])
                for wk in added:
                    if wk not in unset:
                        set_to_process.append(wk)
                        
            tss = {}
            for q in J.gens():
                qm = self.lm(q)
                if qm.is_constant():
                    continue
                for vk in vkstart:
                    v = mod.get_vss(vk[0])[vk[1]][vk[2]]
                    if self.R0.monomial_divides(liftR(qm,self.R0),v.lm()):
                        tss.setdefault(vk,[]).append(q.coefficient(qm))
                        
            vk = min(vkstart,key=lambda vk: (sum(vk[0]), max(vk[0]),
                                             self.tm_dim_tperp(vk[0])-num_in[vk[0]]-num_out[vk[0]],
                                             params_needed.get(vk,0),
                                             len(tss.get(vk,[]))))
            
            set_stack.append(vk)
            vkremove(vk)
            height[0] += 1
            try:
                tm,wt,vi = vk
                vs = mod.get_vss(tm)[wt]
                v = vs[vi]
                
                try:
                    var_lp_val = next(v for x,v in bp.extend() if x == mod.get_vss(tm)[wt][vi].lm())
                except StopIteration:
                    var_lp_val = 0.5
                
                assert vk in unset
                assert bp.get_min(v.lm()) == 0
                assert bp.get_max(v.lm()) == 1
                
                def printstep(extra=''):
                    remhi = len(mod.get_vs(tm)) - self.tm_target_out(tm) - num_in[tm]
                    remlo = self.tm_target_out(tm) - num_out[tm]
                    vstr = str(v.lm())
                    vstr_width[0] = max(vstr_width[0],len(vstr))
                    print ('%s%s %s xi=%d %.3f %s%d %d %s' % ('sound ' if sound else '',
                               str(tm),vstr.rjust(vstr_width[0]),xi,var_lp_val,
                                  ' '*height[0],remhi,remlo,extra))
                # printstep('all')
                
                ts = tss.get(vk,[])
                if len(ts) > 0:
                    print (ts)
                    
                def getJs(J,S,ts,leave_in):
                    def dfscur(I,S,ts):
                        if len(ts) == 0:
                            yield (I,S,False)
                            return
                        t = ts[0]
                        if t in I:
                            for IS in dfscur(I,S,ts[1:]):
                                yield IS
                            return
                        for Icur in (I+[t]).minimal_associated_primes():
                            if any(s in Icur for s in S):
                                continue
                            for IS in dfscur(Icur,S,ts[1:]):
                                yield IS
                        yield (I,S+[t],True)
                    if J.ring().ngens() == self.nvarsx:
                        # no parameters
                        yield (J,S,len(ts) > 0)
                        return
                    I,coerce = self.get_param_ideal(J)
                    S = list(map(coerce,S))
                    ts = list(map(coerce,ts))
                    ts.sort(key=lambda t: t.degree())
                    for I,S,curfull in dfscur(I,S,ts) if leave_in else \
                            ((Ic,S,False) for Ic in (I+ts).minimal_associated_primes() 
                             if not any(s in Ic for s in S)):
                        JI = gb(J+[self.lift_param(J.ring(),t) for t in I.gens()])
                        yield (JI, [self.lift_param(J.ring(),s) for s in S], curfull)
                        
                def rec(J,S,xi,curfull):
                    in_lti = guaranteed_lms(J,S)
                    if not sound: 
                        if not curfull and in_lti(v.lm()):
                            # If killing ts causes a more leading term of
                            # another generator to zero out leaving an
                            # intermediate term to lead, which is us. This
                            # doesn't catch the case where it remains
                            # generically present, but still may only over
                            # approximate the resulting ideals
                            return
                        if any(in_lti(w.lm()) for w in vs[:vi] if bp.get_max(w.lm()) == 0):
                            # Check similarly as above for smaller monomials of
                            # this same weight we have previously excluded (common case)
                            return
                    sav = set()
                    for wk in unset:
                        ws = mod.get_vss(wk[0])[wk[1]]
                        p = ws[wk[2]].lm()
                        if in_lti(p):
                            bp.set_min(p,1)
                            lpset.add((p,1))
                            sav.add(wk)
                            num_in[wk[0]] += 1
                            if wk[0] in tms:
                                for j in range(wk[2]+1,len(ws)):
                                    params_needed[(wk[0],wk[1],j)] -= 1
                                    assert params_needed[(wk[0],wk[1],j)] >= 0
                    for wk in sav:
                        unset.remove(wk)
                    if not curfull:
                        bp.set_max(v.lm(),0)
                        lpset.add((v.lm(),0))
                        unset.remove(vk)
                        num_out[tm] += 1
                    assert bp.get_min(v.lm()) == (1 if curfull else 0)
                    assert bp.get_max(v.lm()) == (1 if curfull else 0)
                    if bp.extend() is not None:
                        if vk in tmsvks:
                            printstep()
                        for JSxi in dfs(J,S,xi):
                            yield JSxi
                    if not curfull:
                        bp.set_max(v.lm(),1)
                        lpset.remove((v.lm(),0))
                        unset.add(vk)
                        num_out[tm] -= 1
                    for wk in sav:
                        ws = mod.get_vss(wk[0])[wk[1]]
                        bp.set_min(ws[wk[2]].lm(),0)
                        lpset.remove((ws[wk[2]].lm(),1))
                        unset.add(wk)
                        num_in[wk[0]] -= 1
                        if wk[0] in tms:
                            for j in range(wk[2]+1,len(ws)):
                                params_needed[(wk[0],wk[1],j)] += 1
                                    
                def addcur(J,S,xi):
                    cur_params = J.ring().ngens() - self.nvarsx
                    yi = xi + params_needed[vk]
                    if yi > cur_params:
                        R = self.get_ring(yi)
                        J = R.ideal([liftR(p,R) for p in J.gens()])
                        S = [liftR(s,R) for s in S]
                    R = J.ring()
                    p = liftR(v,R) + R.sum([liftR(w,R)*R.gen(self.nvarsx+xi+wi)
                                for wi,w in enumerate(w for w in mod.get_vss(tm)[wt][:vi] if bp.get_min(w.lm()) == 0)])
                    for JSxi in rec(gb(J+p),S,yi,True):
                        yield JSxi
                
                if var_lp_val > 0.5 and vk in tmsvks and preset_choices.get(v.lm(),True):
                    for JSxi in addcur(J,S,xi):
                        yield JSxi
                        
                for Jt,St,curfull in getJs(J,S,ts,leave_in=vk not in tmsvks):
                    if preset_choices.get(v.lm(),curfull) == curfull:
                        for JSxi in rec(Jt,St,xi,curfull):
                            yield JSxi
                        
                if var_lp_val <= 0.5 and vk in tmsvks and preset_choices.get(v.lm(),True):
                    for JSxi in addcur(J,S,xi):
                        yield JSxi
                        
            finally:
                height[0] -= 1
                while len(set_stack) > 0:
                    vk = set_stack.pop()
                    vkadd(vk)
                        
        num = 0
        for J,S,xi in dfs(J,S,xi):
            if precise_gb_bound:
                J = J.ring().ideal([p for p in J.gens() if not p.is_zero() and 
                    all([e<=f for e,f in zip(self.polynomial_tm(p),tmmax)])])
                J += [liftR(p,J.ring()) for p in Jorig.gens()]
            if xi == 0 or sound:
                yield self.quick_simplify(J,S,xi)
            else:
                # yield (J,S,xi)
                # continue
                for JSxi in self.get_ideals_for_tm(tms,tmcs,J,S,xi,search_above=search_above,
                               lp_interval_length_bound=lp_interval_length_bound,precise_lp=precise_lp,mod=mod,sound=True):
                    yield JSxi
            num += 1
        if num == 0:
            print ("get_ideals_for_tm: no solution found")
            # ttt.append((J,S,xi,tms,1))
#            if len(tms) > 0:
#                write_JS(J,S,'bapolarsols/mmult2nosol/nosol%d_%s.txt'% (cnt[0],str(tms[0])))
#                cnt[0] += 1

    def tm_dim(self, tm):
        return prod(binomial(d+t-1,t) for t,d in zip(tm,self.dims))

    def tm_dim_tperp(self, tm):
        if all(a <= b for a,b in zip(tm,self.Tdeg)):
            return len(self.zero.get_vs(tm))
        else:
            return self.tm_dim(tm)

    def tm_codim_tperp(self, tm):
        return self.tm_dim(tm) - self.tm_dim_tperp(tm)
    
    def tm_target_out(self, tm): 
        return min(ba.tm_dim_tperp(tm),ba.r-ba.tm_codim_tperp(tm))

    def polynomial_tm(self, p):
        if p.is_zero():
            return (0,)*len(self.Tdeg)
        e = p.exponents()[0]
        return tuple(sum(k for _,k in e[a:b].sparse_iter()) 
                for a,b in zip(self.vixs,self.vixs[1:]))

    def polynomial_weight(self, p):
        e = p.exponents()[0]
        if e.is_constant():
            return (0,)*self.vwts.nrows()
        return tuple(sum(k*self.vwts[:,xi] for xi,k in e.sparse_iter()).list())

    def all_monomials(self, tm, R=None):
        if R is None:
            R = self.get_ring()
        for y in cartesian_product([IntegerVectors(a,d) 
                for a,d in zip(tm,self.dims) ]):
            yield R.prod(x**k for k,x in 
                    zip((k for e in y for k in e),R.gens()) if k>0)

    def bounded_gb(self, J, d):
        from sage.libs.singular.option import opt,opt_ctx
        twostd = sage.libs.singular.function_factory.ff.twostd
        sav = [p for p in J.gens() if p.degree() > d]
        with opt_ctx:
            opt['deg_bound'] = d
            return J.ring().ideal(twostd(J)+sav,side='twosided')
    
    # For testing purposes
    # This gives the monomials not in the leading term ideal for a generic point
    # in the family of ideals corresponding to J. If the family J is not
    # irreducible, this gives the intersection of this set of monomials over all
    # the irreducible components of J
    # it is assumed that the family J is nonempty (1 not in J)
    def generic_codim(self,J,tm):
        msin = set()
        for p in J.gens():
            m = self.lm(p,True)
            if not m.is_constant(): # otherwise a defining equation of the family
                for n in self.all_monomials(map(sub,tm,self.polynomial_tm(m))):
                    msin.add(m*n)
        return [m for m in self.all_monomials(tm) if m not in msin]

    def simple_greedy_extend(self, J, tm):
        assert J.ring() is self.get_ring()
        cod = len(self.generic_codim(J,tm))
        assert cod >= self.r
        for v in ba.get_vorder(tm)[0]:
            if v not in J:
                K = (J + v).twostd()

                codcur = len(self.generic_codim(K,tm))
                if codcur < self.r:
                    continue
                ok = True
                for s in range(len(tm)):
                    if len(self.generic_codim(K, tm[:s]+(tm[s]+1,)+tm[s+1:])) < self.r:
                        ok = False
                        break
                if not ok:
                    continue
                J = K
                if codcur == self.r:
                    return J

    def sum_ideals(self,Jxis):
        sxis = sum(xi for J,xi in Jxis)
        R = self.get_ring(sxis)
        gens = []
        xcum = 0
        for J,xi in Jxis:
            gens.extend([p(R.gens()[:self.nvarsx] +
                R.gens()[self.nvarsx+xcum:self.nvarsx+xcum+xi] +
                (R.zero(),) * (J.ring().ngens()-self.nvarsx-xi))
                for p in J.gens()])
            xcum += xi
        return (R.ideal(gens,side='twosided'),xcum)
    
    def sum_ideals_S(self,Jxis):
        sxis = sum(xi for J,S,xi in Jxis)
        R = self.get_ring(sxis)
        gens = []
        Sgens = []
        xcum = 0
        for J,S,xi in Jxis:
            gens.extend([p(R.gens()[:self.nvarsx] +
                R.gens()[self.nvarsx+xcum:self.nvarsx+xcum+xi] +
                (R.zero(),) * (J.ring().ngens()-self.nvarsx-xi))
                for p in J.gens()])
            Sgens.extend([ p(R.gens()[:self.nvarsx] +
                R.gens()[self.nvarsx+xcum:self.nvarsx+xcum+xi] +
                (R.zero(),) * (J.ring().ngens()-self.nvarsx-xi))
                for p in S ])
            xcum += xi
        return (R.ideal(gens,side='twosided'),Sgens,xcum)

    def quick_simplify(self,J,S=None,xi=None,simple=False):
        if xi is None:
            xi = J.ring().ngens() - self.nvarsx
        if xi == 0:
            return (J,[J.ring().one()],0)
        
        # if variables never occur, add a relation for them so they will be
        # eliminated in the next step. Missing variables should only be
        # produced with sound=True searches
        tis = set()
        for p in J.gens():
            for e in p.exponents(): 
                for ti,_ in e[self.nvarsx:].sparse_iter():
                    tis.add(ti)
        J += [J.ring().gen(self.nvarsx+ti) for ti in set(range(xi)).difference(tis)]
        
        I,coerce = self.get_param_ideal(J)
        subs = list(I.ring().gens()[:xi]) + [I.ring().zero()]*(I.ring().ngens()-xi)
        for ti,t in enumerate(I.ring().gens()[:xi]):
            subs[ti] = I.reduce(t)
        if not simple and xi > 1: # try to write each parameter in terms of the rest
            I = ideal([p(subs) for p in I.gens()])
            for ti,t in enumerate(I.ring().gens()[:xi]):
                if subs[ti] == t:
                    R = PolynomialRing(I.base_ring(),(t,)+tuple(I.ring().gens()[:ti]+I.ring().gens()[ti+1:]),
                               order = TermOrder('degrevlex',1)+ TermOrder('degrevlex',len(subs)-1))
                    tR = R(t)
                    tr = I.change_ring(R).reduce(tR)
                    if tr != tR:
                        subs[ti] = I.ring()(tr)
                        I = ideal([p(subs) for p in I.gens()])
        tis = set()
        for p in subs:
            for e in p.exponents():
                for ti,_ in e.sparse_iter():
                    tis.add(ti)
        tis = sorted(tis)
        ssubs = [I.ring().zero()]*I.ring().ngens()
        for j,ti in enumerate(tis):
            ssubs[ti] = I.ring().gen(j)
        subs = [t(*ssubs) for t in subs]

        R = self.get_ring(len(tis))
        if S is None:
            S = [R.one()]
        else:
            S = [self.lift_param(R,I.reduce(coerce(s)(subs))) for s in S]
            S = [R.one()] + [s for s in S if not s.is_constant()]
        
        subs = list(R.gens()[:self.nvarsx]) + [self.lift_param(R,t) for t in subs]
        Jg = [p(*subs) for p in J.gens()]
        Jg = [p for p in Jg if not p.is_zero()]
        J = R.ideal(Jg,side='twosided')
        return (J,S,len(tis))

    def leading_term_ideal(self,J,S):
        R = PolynomialRing(self.F,self.R0.gens()[:self.nvars])
        ltIlb = []
        ltIub = []
        for p in J.gens():
            m = self.lm(p)
            if not m.is_constant():
                t = p.coefficient(m)
                m = m(R.gens()+(R.zero(),)*(p.parent().ngens()-R.ngens()))
                ltIub.append(m)
                Jt = (J+t).twostd()
                if any(s in Jt for s in S):
                    ltIlb.append(m)
        return R.ideal(ltIlb), R.ideal(ltIub)
    
    def ideal_constant_generators(self,J):
        return self.R0.ideal([liftR(p,self.R0) for p in J.gens() if 
              all(e[self.nvars:].is_constant() for e in p.exponents())],side='twosided')
                        
    def target_hilbert_series(self):
        ts = FractionField(PolynomialRing(QQ,'t',len(self.dims))).gens()
        hs = self.r/prod(1-t for t in ts)
        k = 0
        while True:
            done = True
            for tm in IntegerVectors(k,len(self.dims)):
                tm = tuple(tm)
                if self.tm_dim(tm) < self.r:
                    hs -= (self.r-self.tm_dim(tm))*prod(t**e for t,e in zip(ts,tm))
                    done = False
            if done:
                break
            k += 1
        return hs

    def target_hilbert_series_td(self):
        t = FractionField(PolynomialRing(QQ,'t')).gen()
        hs = self.target_hilbert_series()
        return hs((t,)*hs.parent().ngens())

    def Tperp(self):
        tms = [(0,)*i+(e+1,)+(0,)*(len(self.Tdeg)-i-1)
             for i,e in enumerate(self.Tdeg)]
        tms.extend(product(*[range(e+1) for e in self.Tdeg]))
        return self.R0.ideal([v for tm in tms for vs in 
                        self.get_vss(tm).values() for v in vs]).twostd()
class JDat:
    def __init__(self, ba, J):
        assert all(e[ba.nvars:].is_constant() for p in J.gens() for e in p.exponents())
        self.ba = ba
        self.R = PolynomialRing(ba.F,ba.R0.gens()[:ba.nvars])
        self.J = self.R.ideal([liftR(p,self.R) for p in J.gens()])
        self.ltJ = self.R.ideal([p.lm() for p in self.J.groebner_basis()])
        self.kbase = cache(self.kbase)
        self.get_vss = cache(self.get_vss)
        self.get_vs = cache(self.get_vs)
        self.get_poset = cache(self.get_poset)
        self.get_poset_info_for_search = cache(self.get_poset_info_for_search)
        self.get_quiver_representation = cache(self.get_quiver_representation)
        self.get_linear_program = cache(self.get_linear_program)
        
    def kbase(self,d):
        return [liftR(m,ba.R0) for m in ff.kbase(self.J,d)]
    
    def get_vss(self,tm):
        vs = self.kbase(sum(tm))
        vss = {}
        for v in vs:
            if ba.polynomial_tm(v) == tm:
                vss.setdefault(ba.polynomial_weight(v),[]).append(v)
        todel = []
        for wt in vss:
            vs = vss[wt]
            vs.sort(reverse=True)
            if all(a <= b for a,b in zip(tm,ba.Tdeg)): # restrict to Tperp
                ms = []
                msix = {}
                eqs = {}
                for vi,v in enumerate(vs):
                    Td = ba.T
                    for i,k in v.exponents()[0].sparse_iter():
                        for _ in range(k):
                            Td = Td.derivative(Td.parent().gen(i))
                    for m,c in Td.dict().items():
                        if m not in msix:
                            msix[m] = len(ms)
                            ms.append(m)
                        eqs[(msix[m], vi)] = c
                eqs = matrix(ba.F, len(ms), len(vs), eqs)
                K = eqs.right_kernel_matrix()
                vs = [sum(e*vs[i] for i,e in enumerate(r)) for r in K]
            vss[wt] = list(reversed(vs))
            if len(vss[wt]) == 0:
                todel.append(wt)
        for wt in todel:
            del vss[wt]
        return vss
    
    def get_vs(self, tm):
        vss = self.get_vss(tm)
        P = self.get_poset((),(tm,))
        return [vss[wt][i] for _,wt,i in P.linear_extension()]
    
    def get_poset(self, tms, tms_high_to_low):
        vsix = {}
        lmtov = {}
        tms = tms + tms_high_to_low
        g = DiGraph()
        for tm in tms:
            for wt,vs in self.get_vss(tm).items():
                for vi,v in enumerate(vs):
                    lmtov[v.lm()] = v
                    vk = (tm,wt,vi)
                    vsix[v] = vk
                    g.add_vertex(vk)
        for tm in tms:
            for vs in self.get_vss(tm).values():
                for x in ba.R0.gens()[ba.nvars:ba.nvars+len(ba.xs)]:
                    wi = 0
                    for v in vs:
                        vk = vsix[v]
                        w = liftR(self.J.reduce(liftR(x*v-v*x,self.R)),ba.R0)
                        if w.is_zero():
                            continue
                        _,wt,wicur = vsix[lmtov[w.lm()]]
                        wi = max(wi,wicur)
                        for i in range(wi+1):
                            g.add_edge(vk,(tm,wt,i))
                        wilast = wi
        for tm in tms_high_to_low:
            for vs in self.get_vss(tm).values():
                for v,w in zip(vs,vs[1:]):
                    g.add_edge(vsix[w],vsix[v])
                            
        tmposet = Poset((tms,lambda a,b: all(i<=j for i,j in zip(a,b))))
        for tm,tmr in tmposet.cover_relations():
            for vs in self.get_vss(tm).values():
                for v in vs:
                    vk = vsix[v]
                    for m in ba.all_monomials(tuple(map(sub,tmr,tm))):
                        w = v*m
                        w = liftR(self.J.reduce(liftR(v*m,self.R)),ba.R0)
                        if w.is_zero():
                            continue
                        wk = vsix[lmtov[w.lm()]]
                        g.add_edge(vk,wk)
                    
        return Poset(g)
    
    def get_poset_info_for_search(self, tms, tms_high_to_low):
        P = self.get_poset(tms,tms_high_to_low)
        wkabove = {wk : [] for wk in P}
        vkbelowcnt = {wk : 0 for wk in P}
        for vk,wk in P.cover_relations():
            wkabove[vk].append(wk)
            vkbelowcnt[wk] += 1
        vkstart = set()
        for vk in P:
            if vkbelowcnt[vk] == 0:
                vkstart.add(vk)
        return vkstart,vkbelowcnt,wkabove
    
    def get_quiver_representation(self,tms,tmcs=()):
        # ps in span of qs, find matrix expressing each 
        # p as a linear combination of the qs
        def get_matrix(qs,ps):
            ms = []
            mdict = {}
            A = {}
            for pi,p in enumerate(ps):
                for m in p.monomials():
                    if m not in mdict:
                        mdict[m] = len(ms)
                        ms.append(m)
                    A[(mdict[m],pi)] = p.monomial_coefficient(m)
            B = {}
            for pi,p in enumerate(qs):
                for m in p.monomials():
                    if m not in mdict:
                        mdict[m] = len(ms)
                        ms.append(m)
                    B[(mdict[m],pi)] = p.monomial_coefficient(m)
            A = matrix(ba.F,len(ms),len(ps),A)
            B = matrix(ba.F,len(ms),len(qs),B)
            return B.solve_right(A)
            
        P = Poset((tms+tmcs,lambda a,b: all(i<=j for i,j in zip(a,b))))
        Q = DiGraph()
        for tm,tmr in P.cover_relations():
            for wt,vs in self.get_vss(tm).items():
                for m in ba.all_monomials(tuple(map(sub,tmr,tm))):
                    vms = [liftR(self.J.reduce(liftR(v*m,self.R)),ba.R0)
                            for v in vs]
                    if all(vm.is_zero() for vm in vms):
                        continue
                    wtup = ba.polynomial_weight(next(vm for vm in vms if not vm.is_zero()))
                    ws = self.get_vss(tmr)[wtup]
                    Q.add_edge((tm,wt),(tmr,wtup),get_matrix(ws,vms))
        for tm in tms:
            for wt,vs in self.get_vss(tm).items():
                for x in ba.R0.gens()[ba.nvars:]:
                    xvs = [liftR(self.J.reduce(liftR(x*v-v*x,self.R)),ba.R0)
                            for v in vs]
                    if all(xv.is_zero() for xv in xvs):
                        continue
                    wtup = ba.polynomial_weight(next(vm for vm in xvs if not vm.is_zero()))
                    Q.add_edge((tm,wt),(tm,wtup),get_matrix(self.get_vss(tm)[wtup],xvs))
        return Q
    
    def get_linear_program(self, tms, tmcs=(), interval_length_bound=6, precise=False):
        lp = MixedIntegerLinearProgram(solver='GLPK')
        # lp.set_binary(lp.default_variable())
        # # if set_real, be sure to change exact comparisons to within machine precision
        lp.set_real(lp.default_variable())
        lp.set_min(lp.default_variable(),0)
        lp.set_max(lp.default_variable(),1)

        for tm in tms:
            lp.add_constraint(lp.sum(1-lp[v.lm()] 
                    for vs in self.get_vss(tm).values() for v in vs) == self.ba.tm_target_out(tm))
        for tm in tmcs:
            lp.add_constraint(lp.sum(1-lp[v.lm()] 
                    for vs in self.get_vss(tm).values() for v in vs) >= self.ba.tm_target_out(tm))

        constraints = SetSystem()
        nconstraints = 0
        def add_constraint(rowvks,colvks,kerdim,deduplp=None):
            assert kerdim <= len(colvks)
            if len(colvks) == kerdim:
                return False
            S = [('r',vk) for vk in rowvks]+[('c',vk) for vk in colvks]
            for S2,(ci,kerdim1) in constraints.iter_sets([],S):
                rowvks1 = []
                colvks1 = []
                for tag,vk in S2:
                    if tag == 'r':
                        rowvks1.append(vk)
                    else:
                        colvks1.append(vk)
                if kerdim-len(colvks) >= kerdim1-len(colvks1):
                    # this condition is one we have added plus conditions 
                    # x <= 1 or 0 <= x
                    return False
            if deduplp is not None:
                from sage.numerical.mip import MIPSolverException
                cc = deduplp.sum(deduplp[self.get_vss(tm)[vk][i].lm()] for tm,vk,i in rowvks) + kerdim - \
                     deduplp.sum(deduplp[self.get_vss(tm)[vk][i].lm()] for tm,vk,i in colvks)
                deduplp.set_objective(cc)
                deduplp.set_min(deduplp.default_variable(),0)
                deduplp.set_max(deduplp.default_variable(),1)
                s = deduplp.solve()
                deduplp.set_objective(0)
                if s >= -1e-9:
                    return False
                deduplp.add_constraint(cc >= 0)
            prows = [self.get_vss(tm)[vk][i] for tm,vk,i in rowvks]
            pcols = [self.get_vss(tm)[vk][i] for tm,vk,i in colvks]
            if len(prows) > 1:
                print(prows,pcols,kerdim)
            nonlocal nconstraints
            constraints.add_set(S, (nconstraints,kerdim))
            nconstraints += 1
            lp.add_constraint(lp.sum(lp[self.get_vss(tm)[vk][i].lm()]
                                         for tm,vk,i in colvks)
                             <= lp.sum(lp[self.get_vss(tm)[vk][i].lm()] 
                                         for tm,vk,i in rowvks) + kerdim)
            return True
            
        Q = self.get_quiver_representation(tms,tmcs)
        def add_relations(to,froms,deduplp,interval_length_bound=oo):
            if len(froms) == 0:
                return
            tosize = len(self.get_vss(to[0])[to[1]])
            sol_ss = {}
            # for l in range(1,tosize+1):
            for l in list(range(1,min(interval_length_bound,tosize-1)+1))+[tosize]:
            # for l in range(1,min(interval_length_bound,tosize)+1):
                for a in range(tosize-l+1):
                    b = a+l
                    sol_ss[(a,b)] = SetSystem()
                    ss = sol_ss[(a,b)]
                    rowvks = [to+(i,) for i in range(a,b)]
                    Ms = []
                    for f in froms:
                        M = Q.edge_label(f,to)
                        # want largest fb so that the subspace indexed up to fb maps into 
                        # the subspace indexed up to b, ie, so that M[b:,:fb].is_zero()
                        fb = next(fb for fb in range(M.ncols(),-1,-1) if M[b:,:fb].is_zero())
                        Ms.append(M[a:b,:fb])
                    r = block_matrix([Ms]).rank()
                    if r != Ms[0].nrows():
                        # could actually continue here with kerdim equal to defect if want
                        continue
                    cols_ss = SetSystem()
                    cols = []
                    cols_ech = []
                    cols_ech_pivots = []
                    jxs = [0 for _ in froms]
                    def dfs():
                        if len(cols) == l:
                            assert l > 0
                            colvks = [froms[ii]+(jx,) for ii,jx in cols]
                            ss.add_set(colvks,None)
                            have = False
                            for l1 in range(1,l):
                                if (a,a+l1) in sol_ss and (a+l1,b) in sol_ss:
                                    ss1 = sol_ss[(a,a+l1)]
                                    ss2 = sol_ss[(a+l1,b)]
                                    if any(True for psol,_ in ss1.iter_sets(Shi=colvks) 
                                           for _ in ss2.iter_sets(Shi=[vk for vk in colvks if vk not in psol])):
                                        have = True
                                        break
                            if not have:
                                add_constraint(rowvks,colvks,0,deduplp)
                            return
                        for ii,jx in enumerate(jxs):
                            iiprev = (len(jxs)+ii-1) % len(jxs)
                            if l > interval_length_bound and any(j != 0 for j in jxs) and \
                                    jx == 0 and jxs[iiprev] < Ms[iiprev].ncols():
                                continue
                            startjx = jx
                            while jx < Ms[ii].ncols():
                                nextcol = Ms[ii][:,jx]
                                for j,col in zip(cols_ech_pivots,cols_ech):
                                    nextcol -= nextcol[j,0]*col
                                assert all(nextcol[j,0] == 0 for j in cols_ech_pivots)
                                if not nextcol.is_zero():
                                    break
                                jx += 1
                            if jx == Ms[ii].ncols():
                                continue
                            nextpivot = next(ei for ei,e in enumerate(nextcol.list()) if e != 0)
                            nextcol /= nextcol[nextpivot,0]
                            cols.append((ii,jx))
                            if any(True for _ in cols_ss.iter_sets([],cols)):
                                cols.pop()
                                continue
                            cols_ech.append(nextcol)
                            cols_ech_pivots.append(nextpivot)
                            jxs[ii] = jx+1
                            dfs()
                            jxs[ii] = startjx
                            cols_ech_pivots.pop()
                            cols_ech.pop()
                            cols_ss.add_set(cols,None)
                            cols.pop()
                    dfs()
            
        for (tmr,wtr) in Q:
            froms = [(tm,wt) for (tm,wt),_,_ in Q.incoming_edges((tmr,wtr))]
            # deduplp = None
            deduplp = MixedIntegerLinearProgram(solver='GLPK',maximization=False)
            add_relations((tmr,wtr), froms, deduplp, interval_length_bound=interval_length_bound)
        print(tms,tmcs,lp.number_of_constraints(),'constraints')
        if precise:
            ms = dict([(next(iter(v.dict().keys())),m) for m,v in lp.default_variable().items()])
            msix = dict([(m,next(iter(v.dict().keys()))) for m,v in lp.default_variable().items()])
            for tm in tms:
                xs = [p.lm() for p in self.get_vs(tm)]
                if len(xs) == 0 or len(xs) == self.ba.tm_dim(tm):
                    continue
                curxis = set(msix[m.lm()] for tmcur in (tm,)+tmcs for m in self.get_vs(tmcur))
                lpcur = MixedIntegerLinearProgram(solver='GLPK')
                for lo,(xis,cs),hi in lp.constraints():
                    if all(xi in curxis for xi in xis):
                        lpcur.add_constraint(lpcur.sum(c*lpcur[ms[xi]] for xi,c in zip(xis,cs)),min=lo,max=hi)
                lpcur.set_real(lpcur.default_variable())
                lpcur.set_min(lpcur.default_variable(),0)
                lpcur.set_max(lpcur.default_variable(),1)
                xs = [p.lm() for p in self.get_vs(tm)]
                from sage.numerical.mip import MIPSolverException

                contains_at_least = SetSystem()
                for size in range(1,3):
                    for ys in combinations(xs,size):
                        print(ys)
                        lpcur.set_objective( lpcur.sum(-lpcur[x] for x in ys) )
                        at_least = -floor(lpcur.solve() + 1e-10)
                        if at_least > 0 and all(sum(at_least2) < at_least for _,at_least2 in contains_at_least.iter_unions(Shi=ys)):
                            print('adding', ys, at_least)
                            lpcur.add_constraint(lpcur.sum(lpcur[x] for x in ys) >= at_least)
                            lp.add_constraint(lp.sum(lp[x] for x in ys) >= at_least)
                            contains_at_least.add_set(ys,at_least)
                        
                
                # sols = SetSystem()
                # def extends_to_solution(psol):
                #     if any(True for _ in sols.iter_sets(Slo=psol)):
                #         return True
                #     for x in psol:
                #         lpcur.set_max(lpcur[x],0)
                #     res = list(islice(lp_integer_points(lpcur,xs,fullsol=False),1))
                #     for x in psol:
                #         lpcur.set_max(lpcur[x],1)
                #     if len(res) == 0:
                #         return False
                #     else:
                #         sol = [x for x,v in res[0].items() if v==0]
                #         sols.add_set(sol)
                #         return True
                # constraints = SetSystem()
                # ix = []
                # def dfs(size_target):
                #     if len(ix) > size_target or any(True for _ in constraints.iter_sets(Shi=ix)):
                #         return
                #     print(size_target,ix)
                #     if extends_to_solution([xs[i] for i in ix]):
                #         for i in range(0 if len(ix) == 0 else ix[-1]+1, len(xs)):
                #             ix.append(i)
                #             dfs(size_target)
                #             ix.pop()
                #     elif len(ix) == size_target:
                #         print('adding',ix)
                #         constraints.add_set(ix)
                #         lp.add_constraint(lp.sum(lp[xs[i]] for i in ix) >= 1)
                #         lpcur.add_constraint(lpcur.sum(lpcur[xs[i]] for i in ix) >= 1)
                # for size_target in range(1, min(2,self.ba.tm_target_out(tm))+1):
                #     dfs(size_target)
                

                # st = []
                # cur = []
                # bads = SetSystem()
                # def dfs(check_solvable=True):
                #     try:
                #         if check_solvable:
                #             lpcur.solve()
                #     except MIPSolverException:
                #         return False
                #     csol = lpcur.get_values(lpcur.default_variable())
                #     remxs = [x for x in xs if lpcur.get_min(lpcur[x]) < lpcur.get_max(lpcur[x])]
                #     if len(remxs) == 0:
                #         return True
                #     x = min(remxs, key=lambda x: abs(csol[x]-round(csol[x])))
                #     print('%s%s' % (''.join(st),str(x)))
                #     # hi_first = csol[x] < 0.5
                #     hi_first = csol[x] >= 0.5
                #     if hi_first:
                #         st.append(' ')
                #         lpcur.set_min(lpcur[x],1)
                #         cur.append((x,True))
                #         solhi = dfs(csol[x] < 1-1e-10)
                #         cur.pop()
                #         lpcur.set_min(lpcur[x],0)
                #         st[-1] = '.'
                #     else:
                #         st.append(' ')
                #     lpcur.set_max(lpcur[x],0)
                #     cur.append((x,False))
                #     sollo = dfs(csol[x] > 1e-10)
                #     cur.pop()
                #     lpcur.set_max(lpcur[x],1)
                #     if not hi_first:
                #         st[-1] = '.'
                #         lpcur.set_min(lpcur[x],1)
                #         cur.append((x,True))
                #         solhi = dfs(csol[x] < 1-1e-10)
                #         cur.pop()
                #         lpcur.set_min(lpcur[x],0)
                #     st.pop()
                #     if sollo != solhi:
                #         bad = cur + [(x,sollo)]
                #         lpcur.set_min(lpcur[x],1 if sollo else 0)
                #         lpcur.set_max(lpcur[x],1 if sollo else 0)
                #         for xcur in lpcur.default_variable().keys():
                #             lpcur.set_min(lpcur[xcur],0)
                #             lpcur.set_max(lpcur[xcur],1)
                #         for xcur,v in bad:
                #             lpcur.set_min(lpcur[xcur],v)
                #             lpcur.set_max(lpcur[xcur],v)
                #         i = 0
                #         while i < len(bad):
                #             xcur,v = bad[i]
                #             lpcur.set_min(lpcur[xcur],0)
                #             lpcur.set_max(lpcur[xcur],1)
                #             # if any(True for _ in lp_integer_points(lpcur,xs)):
                #             #     lpcur.set_min(lpcur[xcur],1 if v else 0)
                #             #     lpcur.set_max(lpcur[xcur],1 if v else 0)
                #             #     i += 1
                #             # else:
                #             #     assert len(bad) == 1 or i < len(bad) - 1
                #             #     if len(bad) == 1:
                #             #         i += 1
                #             #     else:
                #             #         bad = bad[:i]+bad[i+1:]
                #             try:
                #                 lpcur.solve()
                #                 lpcur.set_min(lpcur[xcur],1 if v else 0)
                #                 lpcur.set_max(lpcur[xcur],1 if v else 0)
                #                 i += 1
                #             except MIPSolverException:
                #                 try:
                #                     assert len(bad) == 1 or i < len(bad) - 1
                #                 except:
                #                     embed()
                #                     raise
                #                 if len(bad) == 1:
                #                     i += 1
                #                 else:
                #                     bad = bad[:i]+bad[i+1:]
                #         if not any(True for _ in bads.iter_sets(Shi=bad)):
                #             print(bad)
                #             bads.add_set(bad)
                #             lpcur.add_constraint(lpcur.sum(lpcur[xcur] if v else 1-lpcur[xcur] for xcur,v in bad) <= len(bad)-1)
                #             lp.add_constraint(lp.sum(lp[xcur] if v else 1-lp[xcur] for xcur,v in bad) <= len(bad)-1)
                #         lpcur.set_min(lpcur[x],0)
                #         lpcur.set_max(lpcur[x],1)
                #         for xcur,v in cur:
                #             lpcur.set_min(lpcur[xcur],1 if v else 0)
                #             lpcur.set_max(lpcur[xcur],1 if v else 0)
                #     return sollo or solhi
                # dfs()
        
        
        #     ms = dict([(next(iter(v.dict().keys())),m) for m,v in lp.default_variable().items()])
        #     msix = dict([(m,next(iter(v.dict().keys()))) for m,v in lp.default_variable().items()])
        #     for tm in tms:
        #         xs = [p.lm() for p in self.get_vs(tm)]
        #         if len(xs) == 0 or len(xs) == self.ba.tm_dim(tm):
        #             continue
        #         curxis = set(msix[m.lm()] for tmcur in (tm,)+tmcs for m in self.get_vs(tmcur))
                
        #         # import islpy as isl
        #         # space = isl.Space.create_from_names(isl.DEFAULT_CONTEXT, set=['x%d' % i for i in curxis])
        #         # constrs = []
        #         # for xi in curxis:
        #         #     constrs.append(isl.Constraint.ineq_from_names(space, {1:1, 'x%d' % xi: -1}))
        #         #     constrs.append(isl.Constraint.ineq_from_names(space, {'x%d' % xi: 1}))
        #         # for lo,(xis,cs),hi in lp.constraints():
        #         #     if all(xi in curxis for xi in xis):
        #         #         if lo is not None:
        #         #             constr = {1: -int(lo)} | {'x%d' % xi : int(c) for xi,c in zip(xis,cs) }
        #         #             constrs.append(isl.Constraint.ineq_from_names(space,constr))
        #         #         if hi is not None:
        #         #             constr = {1: int(hi)} | {'x%d' % xi : -int(c) for xi,c in zip(xis,cs) }
        #         #             constrs.append(isl.Constraint.ineq_from_names(space,constr))
        #         # bset = isl.BasicSet.universe(space).add_constraints(constrs)
        #         # # return bset, 
        #         # tmxis = [msix[m] for m in xs]
        #         # print('projecting...')
        #         # bset = bset.project_out_except(['x%d' % xi for xi in tmxis],[isl.dim_type.out]*len(tmxis))
        #         # print(bset)
        #         # print(tmxis)
        #         # ixs = [s.get_var_dict()['x%d'%xi][1] for xi in tmxis]
        #         # R = PolynomialRing(QQ,'x',len(xs))
        #         # bad_condition = R.ideal(1)
        #         # def process_point(pt):
        #         #     print(pt)
        #         #     nonlocal bad_condition
        #         #     bad_condition = bad_condition.intersection(R.ideal(
        #         #         [R.gen(i) for i in ixs if pt.get_coordinate_val(isl.dim_type.set,i).is_zero()]))
        #         # bset.foreach_point(process_point)
                                        
        #         lp = MixedIntegerLinearProgram(solver='GLPK')
        #         for lo,(xis,cs),hi in lp.constraints():
        #             if all(xi in curxis for xi in xis):
        #                 lp.add_constraint(lp.sum(c*lp[ms[xi]] for xi,c in zip(xis,cs)),min=lo,max=hi)
        #         lp.set_real(lp.default_variable())
        #         lp.set_min(lp.default_variable(),0)
        #         lp.set_max(lp.default_variable(),1)
        #         R = PolynomialRing(QQ,'x',len(xs))
        #         bad_condition = R.ideal(1)
        #         for si,sol in enumerate(lp_integer_points(lp,xs,fullsol=False)):
        #             print(tm,si)
        #             sol = [m for m,v in sol.items() if v==0]
        #             bad_condition = bad_condition.intersection(R.ideal([
        #                 R.gen(xs.index(m)) for m in sol]))
                
        #         print(bad_condition)
        #         for bad in bad_condition.gens():
        #             assert len(bad.exponents()) == 1
        #             assert all(k==1 for i,k in bad.exponents()[0].sparse_iter())
        #             lp.add_constraint(lp.sum(lp[xs[i].lm()] for i,_ in bad.exponents()[0].sparse_iter()) <= bad.degree()-1)
                    
        return lp

def liftR(p,R):
    if p.parent() is R:
        return p
    if p.is_zero():
        return R.zero()
    return p(R.gens()[:p.parent().ngens()]+(R.zero(),)*(p.parent().ngens()-R.ngens()))

def basis_of_lie_algebra(xs):
    n = xs[0].dimensions()[0]
    F = xs[0].base_ring()
    B = copy(xs)
    Q = copy(xs)
    have = matrix(F,len(B),n*n,[x.list() for x in xs],sparse=True)
    have.echelonize()
    assert have.rank() == have.nrows()
    while len(Q) > 0:
        y = Q.pop()
        for x in xs:
            z = x*y-y*x
            havez = block_matrix([[have],[matrix(z.list(),sparse=True)]])
            havez.echelonize()
            if havez.rank() > have.nrows():
                have = havez
                B.append(z)
                Q.append(z)
    return B

def irrelevant_tm_generators(tms):
    R = PolynomialRing(QQ,'x',len(tms[0]))
    Is = [R.ideal([t**(e+1) for e,t in zip(tm,R.gens())]) for tm in tms]
    I = Is[0].intersection(*Is[1:])
    return [tuple(m.exponents()[0]) for m in I.gens()]

def write_JS(J,S,f):
    print ('writing %s' %f)
    f = open(f,'w')
    for s in S:
        f.write(str(s)+'\n')
    f.write('\n')
    for p in J.gens():
        f.write(str(p)+'\n')
    f.close()
    
def read_JS(ba,f):
    f = open(f).read()
    import re
    ti = re.findall('t([0-9]+)',f)
    tmax = 0 if len(ti) == 0 else max(map(eval,ti))+1
    R = ba.get_ring(tmax)
    locals = {str(x):x for x in R.gens()}
    f = f.split('\n')
    S = []
    J = []
    inS = True
    for l in f:
        if len(l) == 0:
            inS = False
            continue
        if inS:
            S.append(R(sage_eval(l,locals=locals)))
        else:
            J.append(R(sage_eval(l,locals=locals)))
    return R.ideal(J,side='twosided'),S,tmax
    
@


\begin{thebibliography}{9}
  \bibitem{levendovsky}
  Levandovskyy, V. Non-commutative Computer Algebra for polynomial algebras:
  Gr\"obner bases, applications and implementation. Doctoral Thesis,
  Universit\"at Kaiserslautern, 2005. Available online at http://kluedo.ub.uni-kl.de/volltexte/2005/1883/.articular 
\end{thebibliography}
\end{document}