def _load_bp(self, fname):
        bp = []
        try:
            with open(fname) as f:
                for line in f:
                    if line.startswith('#'):
                        continue
                    bp_json = json.loads(line)
                    for step in bp_json['steps']:
                        bp.append(
                            Layer(int(step['position']), matrix(step['0']), matrix(step['1'])))

                    assert len(bp_json['outputs'])    == 1 and \
                           len(bp_json['outputs'][0]) == 2
                    first_out = bp_json['outputs'][0][0].lower()
                    if first_out not in ['false', '0']:
                        if first_out not in ['true', '1']:
                            print('warning: interpreting %s as a truthy output' % first_out)
                        bp[-1].zero.swap_columns(0,1)
                        bp[-1].one .swap_columns(0,1)
                    return bp
        except IOError as e:
            print(e)
            sys.exit(1)
        except ValueError as e:
            print('expected numeric position while parsing branching program JSON')
            print(e)
            sys.exit(1)

def run_example(filepath,more_info=False):
    e = E.Example(filepath)
    mA = S.matrix(S.QQ,lmatrix_to_numbers(e.matrix_A))
    mB_s = S.matrix(S.QQ,lmatrix_to_numbers(e.matrix_B_strict))
    mB_w = S.matrix(S.QQ,lmatrix_to_numbers(e.matrix_B_weak))
    print "Checking termination for example '"+ e.example_name +"'"
    print "published in " + e.published_in
    if more_info:
        print "Matrix A:"
        print mA
        print "Matrix jnf(A):"
        print mA.jordan_form(S.QQbar)
        print "Matrix B strict:"
        print mB_s
        print "Matrix B weak:"
        print mB_w
        print ("applicable to complexity theorem: " + str(C.are_absolute_eigenvalues_of_jordan_blocks_distinct(mA)))
        print ("polynomial update:                " + str(C.has_polynomial_growth(mA)))
        print ("max growth:                       " + str(C.pretty_growth(C.max_growth(mA))))
    sys.stdout.flush()
    start_time = time.time()
    result = T.termination_check(mA,mB_s,mB_w,more_info)
    end_time = time.time()
    print "result:",result
    print ("time: %0.4f seconds" % (end_time - start_time))
    print "-"*40+"\n"
    sys.stdout.flush()

    def compatibility_degree(self, alpha, beta):
        if self.is_finite():
            tube_contribution = -1
        elif self.is_affine():
            gck = self.gamma().associated_coroot()
            if any([gck.scalar(alpha) != 0, gck.scalar(beta) != 0]):
                tube_contribution = -1
            else:
                sup_a = self._tube_support(alpha)
                sup_b = self._tube_support(beta)
                if all([x in sup_b for x in sup_a]) or all([x in sup_a for x in sup_b]):
                    tube_contribution = -1
                else:
                    nbh_a = self._tube_nbh(alpha)
                    tube_contribution = len([ x for x in nbh_a if x in sup_b ])
        else:
            raise ValueError("compatibility degree is implemented only for finite and affine types")
        
        initial = self.initial_cluster()
        if alpha in initial:
            return max(beta[initial.index(alpha)],0)

        alphacheck = alpha.associated_coroot()

        if beta in initial:
            return max(alphacheck[initial.index(beta)],0)

        Ap = -matrix(self.rk, map(lambda x: max(x,0), self.b_matrix().list() ) )
        Am =  matrix(self.rk, map(lambda x: min(x,0), self.b_matrix().list() ) )

        a = vector(alphacheck)
        b = vector(beta)

        return max( -a*b-a*Am*b, -a*b-a*Ap*b, tube_contribution )

    def coxeter(self):
        r"""
        Returns a list expressing the coxeter element corresponding to self._B
        (twisted) reflections are applied from top of the list, for example
        [2, 1, 0] correspond to s_2s_1s_0

        Sources == non positive columns == leftmost letters
        """
        zero_vector = vector([0 for x in range(self.rk)])
        coxeter = []
        B = copy(self.B0)
        columns = B.columns()
        source = None
        for j in range(self.rk):
            for i in range(self.rk):
                if all(x <=0 for x in columns[i]) and columns[i] != zero_vector:
                    source = i
                    break
            if source == None:
                if B != matrix(self.rk):
                    raise ValueError("Unable to find a Coxeter element representing self.B0")
                coxeter += [ x for x in range(self.rk) if x not in coxeter]
                break
            coxeter.append(source)
            columns[source] = zero_vector
            B = matrix(columns).transpose()
            B[source] = zero_vector
            columns = B.columns()
            source = None
        return tuple(coxeter)

def kernel_lattice(A, mod=None):
    ''' Lattice of vectors x with Ax = 0 (potentially mod m) '''
    A = matrix(ZZ if mod is None else Integers(mod), A)
    L = [vector(ZZ, row) for row in A.right_kernel().basis()]
    if mod is not None:
        cols = len(L[0])
        for i in range(cols):
            L.append([0]*i + [mod] + [0]*(cols-i-1))
    return matrix(L)

def _blocks_to_quad_form(blcs, p):
    h = matrix([[QQ(0), QQ(1) / QQ(2)],
                [QQ(1) / QQ(2), QQ(0)]])
    y = matrix([[QQ(1), QQ(1) / QQ(2)],
                [QQ(1) / QQ(2), QQ(1)]])
    mat_dict = {"h": h, "y": y}
    mats_w_expt = [(expt, mat_dict[qf] if qf in ("h", "y") else matrix([[qf]]))
                   for expt, qf in blcs]
    qfs = [QuadraticForm(ZZ, m * ZZ(2) * p ** expt) for expt, m in mats_w_expt]
    return reduce(operator.add, qfs)

def _rankin_cohen_triple_det_sym2_pol(k1, k2, k3):
    (r11, r12, r22, s11, s12, s22, t11, t12, t22), (u1, u2) = _triple_gens()

    m0 = matrix([[r11, s11, t11], [2 * r12, 2 * s12, 2 * t12], [k1, k2, k3]])

    m1 = matrix([[r11, s11, t11], [k1, k2, k3], [r22, s22, t22]])

    m2 = matrix([[k1, k2, k3], [2 * r12, 2 * s12, 2 * t12], [r22, s22, t22]])
    Q = m0.det() * u1**2 - 2 * m1.det() * u1 * u2 + m2.det() * u2**2
    return Q

 def matrix_representaion(self, lin_op):
     '''Let lin_op(f, t) be an endomorphsim of self, where f is
     a modular form and t is a object corresponding to a matrix.
     This medthod returns the matrix representation of lin_op.
     '''
     basis = self.basis()
     lin_indep_tuples = self.linearly_indep_tuples()
     m1 = matrix([[f[t] for t in lin_indep_tuples] for f in basis])
     m2 = matrix([[lin_op(f, t) for t in lin_indep_tuples]
                  for f in basis])
     return (m2 * m1 ** (-1)).transpose()

def cvp_embed(L, v, b=None):
    if not b:
        b = max(max(row) for row in L.rows())

    L2 = matrix([list(row) + [0] for row in L] + [list(v) + [b]])
    res = None
    for x in lll(matrix(L2)):
        if x[-1] > 0: x = -x
        if x[-1] == -b:
            u = vector(x[:-1]) + v
            assert in_lattice(L, u)
            if res is None or (v - u).norm() < (v - res).norm():
                res = u
    return res

def mod_right_kernel(A, mod):
    # from https://ask.sagemath.org/question/33890/how-to-find-kernel-of-a-matrix-in-mathbbzn/
    # too slow though
    Zn = ZZ**A.ncols()
    M = Zn/(mod*Zn)
    phi = M.hom([M(a) for a in A] + [M(0) for _ in range(A.ncols()-A.nrows())])
    return matrix([M(b) for b in phi.kernel().gens()])

 def test_division_generators(self):
     prec = 6
     div_consts = [c for c in gens_consts if isinstance(c, ConstDivision)]
     consts = (even_gen_consts() + odd_gen_consts() +
               [CVH(_wt18_consts[0], 2), CVH(sym10_19_consts[0], 2)])
     calculator = CalculatorVectValued(consts, data_dir)
     calculator.calc_forms_and_save(prec, verbose=True, force=True)
     gens_dct = calculator.forms_dict(prec)
     for c in div_consts:
         k = c.weight()
         print "checking when k = %s" % (str(k), )
         if k % 2 == 0:
             sccst = _find_const_of_e4_e6_of_same_wt(18 - k)
             M = Sym10EvenDiv(sccst, prec)
         else:
             sccst = _find_const_of_e4_e6_of_same_wt(19 - k)
             M = Sym10OddDiv(sccst, prec)
         pl = _anihilate_pol(k, M)
         hol_basis = M.basis_of_subsp_annihilated_by(pl)
         N = Sym10GivenWtBase(prec, k, hol_basis)
         # Check this prec is sufficient.
         mt = matrix(QQ, [[b[t] for b in N.basis()]
                          for t in N.linearly_indep_tuples()])
         self.assertTrue(
             mt.is_invertible(), "False when k = %s" % (str(k),))
         # Check our construction gives a holomorphic modular form
         self.assertTrue(N.contains(gens_dct[c]),
                         "False when k = %s" % (str(k),))

    def is_independent(self, v):
        """
        Return True if the Hecke operators in v are independent.

        INPUT:

            - `v` -- four elements of the Hecke algebra mod 2 (represented as matrices)

        OUTPUT:

            - bool

        EXAMPLES::

            sage: from mdsage import *
            sage: C = KamiennyCriterion(29)
            sage: C.is_independent([C.T(1), C.T(2), C.T(3), C.T(4)])
            True
            sage: C.is_independent([C.T(1), C.T(2), C.T(3), C.T(1)+C.T(3)])
            False        
        """
#        X = matrix(GF(2), 4, sum([a.list() for a in v], []))
#        c = sage.matrix.matrix_modn_dense.Matrix_modn_dense(X.parent(),X.list(),False,True)
#        return c.rank() == 4

        # This crashes!  See http://trac.sagemath.org/sage_trac/ticket/8301
        return matrix(GF(2), len(v), sum([a.list() for a in v], [])).rank() == len(v)
        raise NotImplementedError

def _arc(p,q,s,**kwds):
    #rewrite this to use polar_plot and get points to do filled triangles
    from sage.misc.functional import det
    from sage.plot.line import line
    from sage.misc.functional import norm
    from sage.symbolic.all import pi
    from sage.plot.arc import arc

    p,q,s = map( lambda x: vector(x), [p,q,s])

    # to avoid running into division by 0 we set to be colinear vectors that are
    # almost colinear
    if abs(det(matrix([p-s,q-s])))<0.01:
        return line((p,q),**kwds)

    (cx,cy)=var('cx','cy')
    equations=[
            2*cx*(s[0]-p[0])+2*cy*(s[1]-p[1]) == s[0]**2+s[1]**2-p[0]**2-p[1]**2,
            2*cx*(s[0]-q[0])+2*cy*(s[1]-q[1]) == s[0]**2+s[1]**2-q[0]**2-q[1]**2
            ]
    c = vector( [solve( equations, (cx,cy), solution_dict=True )[0][i] for i in [cx,cy]] )
    
    r = norm(p-c)

    a_p, a_q, a_s = map(lambda x: atan2(x[1],x[0]), [p-c,q-c,s-c])
    a_p, a_q = sorted([a_p,a_q])
    if a_s < a_p or a_s > a_q:
        return arc( c, r, angle=a_q, sector=(0,2*pi-a_q+a_p), **kwds)
    return arc( c, r, angle=a_p, sector=(0,a_q-a_p), **kwds)

def isom(A,B):
    # First check that A is a symmetric matrix.
    if not matrix(A).is_symmetric():
        return False
    # Then check A against the viable database candidates.
    else:
        n=len(A[0])
        m=len(B[0])
        Avec=[]
        Bvec=[]
        for i in range(n):
            for j in range(i,n):
                if i==j:
                    Avec+=[A[i][j]]
                else:
                    Avec+=[2*A[i][j]]
        for i in range(m):
            for j in range(i,m):
                if i==j:
                    Bvec+=[B[i][j]]
                else:
                    Bvec+=[2*B[i][j]]
        Aquad=QuadraticForm(ZZ,len(A[0]),Avec)
    # check positive definite
        if Aquad.is_positive_definite():
            Bquad=QuadraticForm(ZZ,len(B[0]),Bvec)
            return Aquad.is_globally_equivalent_to(Bquad)
        else:
            return False

def do_import(ll):
    dim,det,level,gram,density,hermite,minimum,kissing,shortest,aut,theta_series,class_number,genus_reps,name,comments = ll
    mykeys = ['dim','det','level','gram','density','hermite', 'minimum','kissing','shortest','aut','theta_series','class_number','genus_reps','name','comments']
    data = {}
    for j in range(len(mykeys)):
        data[mykeys[j]] = ll[j]
	
    blabel = base_label(data['dim'],data['det'],data['level'], data['class_number'])
    data['base_label'] = blabel
    data['index'] = label_lookup(blabel)
    label= last_label(blabel, data['index'])
    data['label'] = label
 
    lattice = lat.find_one({'label': label})

    if lattice is None:
        print "new lattice"
        print "***********"
        print "check for isometries..."
        A=data['gram'];
        n=len(A[0])
        d=matrix(A).determinant()
        result=[B for B in lat.find({'dim': int(n), 'det' : int(d)}) if isom(A, B['gram'])]
        if len(result)>0:
            print "... the lattice with base label "+ blabel + " is isometric to " + str(result[0]['gram'])
            print "***********"
        else:
            lattice = data
    else:
        print "lattice already in the database"
        lattice.update(data)
    if saving:
        lat.update({'label': label} , {"$set": lattice}, upsert=True)

def sqcap_mul(A, B, n, p, q):
    '''
    Let A and B be square matrices of size
    binomial(n, p) and binomial(n, q).
    Return sqcap multiplication defined in [Bö] as a square matrix of size
    binomial(n, p + q).
    '''
    # if p or q is zero, return immediately.
    if p == 0:
        return B
    elif q == 0:
        return A

    p_dct = _index_dct(n, p)
    q_dct = _index_dct(n, q)
    p_q_dct = _index_dct(n, p + q)

    p_q_lst = permutations_increasing(n, p + q)
    res = matrix([[A.base_ring()(0) for _ in p_q_lst] for _ in p_q_lst])
    for ad in _permutations_increasing(n, p):
        for bd in _permutations_increasing(n, p):
            for add in _permutations_increasing(n, q):
                for bdd in _permutations_increasing(n, q):
                    if all(a not in add for a in ad) and all(b not in bdd for b in bd):
                        a = _concat(ad, add)
                        b = _concat(bd, bdd)
                        s = (_sign(ad, add) * _sign(bd, bdd) *
                             A[p_dct[bd], p_dct[ad]] *
                             B[q_dct[bdd], q_dct[add]])
                        res[p_q_dct[b], p_q_dct[a]] += s
    return binomial(p + q, p) ** (-1) * res

def check_add_qexp(dim, min_det=1, max_det=None, fix=False):
    count = 0
    query = {}
    query['dim'] = int(dim)
    query['det'] = {'$gte' : int(min_det)}
    if max_det:
        query['det']['$lte'] = int(max_det)
    else:
        max_det = "infinity"
    lat_set = lat.find(query)
    print("%s lattices to examine of dimension %s and determinant between %s and %s."
          % (lat_set.count(), dim, min_det, max_det))
    if lat_set.count() == 0:
        return None
    print("checking whether the q expansion is stored...")
    for l in lat_set:
        print("Testing lattice %s" % l['label'])
        if l['theta_series'] == "":
            print("q expansion NOT stored")
            if fix:
                M=l['gram']
                exp=[int(i) for i in gp("Vec(1+2*'x*Ser(qfrep("+str(gp(matrix(M)))+",150,0)))")]
                lat.update({'theta_series': l['theta_series']}, {"$set": {'theta_series': exp}}, upsert=True)
                print("Fixed lattice %s" % l['label'])
        else:
            print("q expansion stored")

def kalmanson_matrix(n, aug=False):
    r,c = triu_indices(n, 1)
    k = len(r)
    inds = np.arange(k)
    row_ind = lambda (i,j): inds[np.logical_and(r==i-1, c==j-1)]
    upright = lambda (i,j): (i,j) if i<j else (j,i)
    get_ind = lambda ind_lst: np.array(map(row_ind, map(upright, grouper(2, ind_lst))))

    rows = []
    for i in range(1,n-2):
        for j in range(i+2, n):
            ind_lst = (i, j+1, i+1, j, i, j, i+1, j+1)
            indices = get_ind(ind_lst)
            row = np.zeros(k, dtype=np.int)
            row[indices] = [-1, -1, 1, 1]           
            rows.append(row)

    sub = len(rows)

    for i in range(2, n-1):
        ind_lst = (i, 1, i+1, n, i, n, i+1, 1)
        indices = get_ind(ind_lst)
        row = np.zeros(k, dtype=np.int)
        row[indices] = [-1, -1, 1, 1]           
        rows.append(row)

    mat = matrix(np.vstack(rows))
    mat.subdivide(sub, None)
    if aug:
        return zero_matrix(len(rows),1).augment(mat)
    else:
        return mat

 def assert_degree_2(self, k):
     es = eisenstein_series_degree2(k, prec=10)
     es1 = sess(weight=k, degree=2)
     self.assertTrue(
         all(es[(n, r, m)] == es1.fourier_coefficient(matrix([[n, ZZ(r) / ZZ(2)],
                                                              [ZZ(r) / ZZ(2), m]]))
             for n, r, m in es.prec))

 def pi(z):
     # Note that f(...) returns a string right now, since it's all done at C level.
     # This will prboably change, breaking this code. 
     M = matrix(k,2,eval(f(z).replace(';',','), {'g':g})).transpose()
     v = M.echelon_form()[0]
     v.set_immutable()
     return v