python-ferari 1.0.0-1 / usr / share / pyshared / ferari / pg.py

# Copyright (C) 2006 Robert C. Kirby
#
# This file is part of FErari.
#
# FErari is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# FErari is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with FErari. If not, see <http://www.gnu.org/licenses/>.
#
# First added:  2005-04-01
# Last changed: 2006-04-01

import numpy,build_tensors,graph
from xpermutations import xuniqueCombinations
from numpy.random import rand
from numpy.linalg import svd, det

precision = 5
eps = .1**precision

def randproj( m , n ):
    """creates a random matrix in R^{m,n}"""
    # maps from R^n to R^m
    if m >= n:
        A = rand( m , n )
        [u,s,vt] = svd( A )
        return u
    else:
        A = rand( n , m )
        [u,s,vt] = svd(A)
        return numpy.transpose(u)

def nnz( u ):
    """Returns the number of items in u that are farther away from
    zero than eps."""
    nnz_cur = 0
    for ui in u:
        if abs( ui ) > eps:
            nnz_cur += 1
    return nnz_cur

def nnz_cmp( x , y ):
    """Allows sorting according to the sparsity of the vectors"""
    nnzx = nnz(x)
    nnzy = nnz(y)
    if nnzx < nnzy:
        return -1
    elif nnzx == nnzy:
        return 0
    else:
        return 1

def compare( u , v ):
    """Lexicographic ordering that respects floating-point
    fuzziness"""
    if len(u) != len(v):
        raise RuntimeError, "can't compare"
    diff = u-v
    for d in diff:
        if abs(d) > eps:
            if d < 0:
                return -1
            else:
                return 1
    return 0

def kvp_cmp( x , y ):
    return compare( x[1] , y[1] )

def normalize_sign( u ):
    """Returns u or -u such that the first nonzero entry (up to eps)
    is positive"""
    for ui in u:
        if abs(ui) > eps:
            if ui < 0.0:
                return -u
            else:
                return u
    return u

def unit_vector( u ):
    mau = max(abs(u))
    if mau < eps:
        raise RuntimeError, "barf: divide by zero"
    return normalize_sign( u / mau )

def filterequal( vs , compare ):
    """vs is a list of key/value pairs such that equal values must be
    adjacent in the list.  Returns a list of lists, each of which is
    the set of key/value pairs corresponding to equal items."""
    outer_index = 0
    equal_items = []
    while outer_index < len( vs ):
        cur_items = [ ]
        inner_index = outer_index
        while inner_index< len(vs) \
                  and compare( vs[inner_index] , \
                               vs[outer_index] ) == 0:
            cur_items.append( vs[inner_index] )
            inner_index += 1
        equal_items.append( cur_items )
        outer_index = inner_index

    return equal_items

def rank( A ):
    """Returns the numerical rank of an array."""
    u,s,vt = svd( A )
    return len( [ si for si in s if abs( si ) > eps ] )

def cohyperplanar(bs,n):
    """Predicate indicating whether the members of bs span
    an n-dimensional space"""
    mat = numpy.array( bs )
    return rank( numpy.array( mat ) ) == n

def strike_col( A , j ):
    m,n = A.shape
    return numpy.take( A , [ k for k in range(0,n) if k != j ] , 1 )

def normalized_cross( vecs ):
    n,d = len( vecs ) , len( vecs[ 0 ] )
    mat = numpy.array( vecs )
    if n != d - 1:
        for v in vecs:
            print v
        raise RuntimeError, "barf"
    # normalize the cross product to have unit value to avoid
    # rounding to zero in the cross product routine.
    foo = numpy.array( \
        [ (-1)**i * det(strike_col(mat,i)) \
          for i in range( d ) ] )
    foo /= max( abs( foo ) )
    return normalize_sign( foo )

def bf_line_finder( vecs , m ):
    tups = [ frozenset( t ) \
             for t in xuniqueCombinations( vecs.keys() , m + 1 ) \
             if cohyperplanar( [ vecs[i] for i in t ] , m ) ]

    gr = dict( [ ( tup , {}) for tup in tups ] )

    for (tup1,tup2) in xuniqueCombinations( tups , 2 ):
        if len( tup1.intersection( tup2 ) ) == m:
            gr[tup1][tup2] = (m,None)
            gr[tup2][tup1] = (m,None)

    ccs = graph.connectedComponents( gr )

    Ls = []
    for cc in ccs:
        Lcur = set()
        for k in cc:
            Lcur.update( k )
        Ls.append( Lcur )

    return [ frozenset( L ) for L in Ls ]

def rp_line_finder( vecs , m ):
    d = len( vecs.itervalues().next() )

    udude = randproj(m+1,d)


    pi_vecs = dict( [ (i,numpy.dot(udude,v ) ) \
                    for (i,v) in vecs.iteritems() ] )

    pi_normals = [ (tup,normalized_cross( [ pi_vecs[t] for t in tup ]))\
                   for tup in xuniqueCombinations( pi_vecs.keys() , \
                                                   m )  ]
    pi_normals.sort( kvp_cmp )

##    for i,n in pi_normals:
##        print i,n

    eitems = filterequal( pi_normals , kvp_cmp )

    Ls = []

    for e in eitems:
        if len( e ) > 1:
            inds = set()
            for (indlist,a) in e:
                inds.update( indlist )
            vs = numpy.array( [ vecs[i] for i in inds ] )
            if rank( vs ) == m:
                Ls.append( frozenset( inds ) )
            else:
                mini_vecs = dict( [ (i,vecs[i]) for i in inds ] )
                Ls.extend( bf_line_finder( mini_vecs , m ) )

    return Ls

def line_graph( Ls , m ):
    """Returns a graph whose elements are the members of Ls and whose
    edges have weights that are the size of the intersection."""
    gr = dict( [ (L,{}) for L in Ls ] )
    for (L1,L2) in xuniqueCombinations(Ls,2):
        n = len( L1.intersection(L2) )
        if n > 0:
            gr[L1][L2] = (m - n,None)
            gr[L2][L1] = (m - n,None)
    return gr

def gen_graph ( vecs , Ls , m ):
    """vecs is a dictionary of label/vector and Ls is a list of sets
    indicating (hyper) lines in a partial geometry based on dependency
    between m+1 things."""
    lg = line_graph( Ls , m )
    gg = {}

    weights = dict( [ (L,m) for L in Ls ] )

    while weights:
        # pick a line with minimal weight
        L = graph.argmin( weights )
        w = weights[L]
        weights.pop( L )

        Ldone = L.intersection( gg )
        # sanity check:
        # if there are k items in Ldone and k <= m
        # then the weight I get had better
        # be m - k
        lldone = len( Ldone )
        if (lldone <= m and w != m - lldone) \
               or (lldone > m and w != 0):
            print ct
            print L
            print Ldone
            print w
            raise RuntimeError, "barf"

        Ltodo = L.difference( Ldone )

        roots = set()
        for r in Ldone:
            if len( roots ) == m:
                break
            else:
                roots.add( r )
        # pick the sparsest new elements as new generator members
	Ltodo_list = list( Ltodo )
	Ltodo_list.sort(nnz_cmp)
        for r in Ltodo_list:
            if len( roots ) == m:
                break
            else:
                roots.add( r )

        done_roots = roots.intersection( gg )
        new_roots = roots.intersection( Ltodo )

        # sanity check: length of new_roots
        # ought to be equal to w
        if len( new_roots ) != w:
            raise RuntimeError, "barf"

        # add new items to gg
        for u in new_roots:
            if u in gg:
                raise RuntimeError, "barf"
            gg[u] = {}

        # new roots go into graph have no out-edges
        # rest of Ltodo go into graph with out-edges to
        # each member of the roots.
        for u in Ltodo.difference( new_roots ):
            if u in gg:
                raise RuntimeError, "barf"
            gg[u] = {}
            for v in roots:
                gg[u][v] = (m,None)

        # update neighboring lines' weight
        for Lnb in lg[L]:
            if Lnb in weights:
                weights[Lnb] = max( 0 , m - len( Lnb.intersection( gg ) ) )

    return gg

def gen_graph_cost( gg , A0dict ):
    """Evaluates the total floating-point cost of the dot product
    algorithm associated with a generation graph."""
    cost = 0
    for u in gg:
        if len( gg[u] ) == 0:
            cost += nnz( A0dict[u] )
        elif len( gg[u] ) == 1:
            cost += gg[u].values()[0][0]
        else:
            cost += len( gg[u] )

    return cost


def process( vecs , mmax ):
    first_index = vecs.iterkeys().next()
    n = len( vecs )
    d = len( vecs[first_index] )
    z = numpy.zeros( (d,) , "d" )

    def filterdict( source , filt ):
        return dict( [ a for a in source.iteritems() \
                       if a[0] not in filt ] )

    parents = {}

    print "finding zeros..."
    zeros = [ x for x in vecs.iteritems() \
              if numpy.alltrue( numpy.allclose( x[1],z,eps ) ) ]

    for z in zeros:
        parents[z[0]] = {}

    remaining = filterdict( vecs , parents )

    print "down to %d vecs" % ( len( remaining ) , )

    print "filtering equal vectors"

    # create the list of projections of vectors, each with positive
    # first nonzero entry
    projected = [ (x[0],normalize_sign(x[1])) \
                  for x in remaining.iteritems() ]

    projected.sort( kvp_cmp )

    eitems = filterequal( projected , kvp_cmp )

    for eit in eitems:
        item_cur = eit[0]
        for item in eit[1:]:
            parents[item[0]] = { item_cur[0] : (0,"e") }

    remaining = filterdict( remaining , parents )

    print "done filtering equal vectors"

    for m in range(1,mmax+1):
        if len( remaining ) <= m:
            break
        print "filtering linear dependence of order " , m
        Ls = rp_line_finder( remaining , m )
        gg = gen_graph( remaining , Ls , m )

        for u in gg:
            if gg[u]:
                parents[u] = dict( [ (v,(m,"lc")) for v in gg[u] ] )

        remaining = filterdict( remaining , parents )

        print "%s vectors remaining" % ( len(remaining,) )

    for r in remaining:
        parents[r] = {}

    return parents

def process2( vecs , max_level ):
    first_index = vecs.iterkeys().next()
    n = len( vecs )
    d = len( vecs[first_index] )
    z = numpy.zeros( (d,) , "d" )

    def filterdict( source , filt ):
        return dict( [ a for a in source.iteritems() \
                       if a[0] not in filt ] )

    parents = {}

    nv = len( vecs )

    zeros = [ x for x in vecs.iteritems() \
              if numpy.alltrue( numpy.allclose( x[1],z,eps ) ) ]

    nz = len( zeros )

    nr = nv - nz

    for z in zeros:
        parents[z[0]] = {}

    remaining = filterdict( vecs , parents )

    projected = [ (x[0],normalize_sign(x[1])) \
                  for x in remaining.iteritems() ]

    projected.sort( kvp_cmp )

    eitems = filterequal( projected , kvp_cmp )

    for eit in eitems:
        item_cur = eit[0]
        for item in eit[1:]:
            parents[item[0]] = { item_cur[0] : (0,"e") }

    remaining = filterdict( remaining , parents )

    ne = nr - len( remaining )
    nr = len( remaining )

    print "done filtering equal vectors"

    nactive = []
    nleft = [nr]
    curcost = []

    for i in range(1,max_level+1):
	print "filtering linear dependence of order " , i
	print "%s vectors remain" % (nleft[-1],)
	Ls = rp_line_finder( remaining , i )

	# figure out who lives in linear relations, who doesn't
	active = set()
	for L in Ls:
	    for u in L:
		active.add( u )

	nactive.append( len( active ) )

	print "%s elements are active" % (len(active),)

	gg = gen_graph( remaining , Ls , i )
	for u in gg:
	    if gg[u]:
		parents[u] = dict(  [(v,(i,"lc")) for v in gg[u] ] )

	remaining = filterdict( remaining , parents )
	nr = len( remaining )
	nleft.append( nr )

	print "%s remain now" % (nr,)

	# temporarily add in everything remaining by brute force to the
	# generation graph, compute the cost, and then remove them

	for r in remaining:
	    parents[r] = {}

	cc = gen_graph_cost( parents , vecs )
	curcost.append( cc )
	print "current cost is down to %s" % (cc,)

	for r in remaining:
	    parents.pop( r )

    for r in remaining:
	parents[r] = {}

    return parents,n,d,nz,ne,nleft,nactive,curcost

def process_gen_graph( vecs ):
	"""Returns the generation graph for vectors active at
	2-dependencies"""
	first_index = vecs.iterkeys().next()
	n = len( vecs )
	d = len( vecs.itervalues().next() )
	z = numpy.zeros( (d,) , "d" )

	def filterdict( source , filt ):
		return dict( [ a for a in source.iteritems() \
                       if a[0] not in filt ] )

	parents = {}

	print "finding zeros..."
	zeros = [ x for x in vecs.iteritems() \
              if numpy.alltrue( numpy.allclose( x[1],z,eps ) ) ]

	for z in zeros:
		parents[z[0]] = {}

	remaining = filterdict( vecs , parents )

	print "down to %d vecs" % ( len( remaining ) , )

	print "filtering equal vectors"

	projected = [ (x[0],normalize_sign(x[1])) \
                  for x in remaining.iteritems() ]

	projected.sort( kvp_cmp )

	eitems = filterequal( projected , kvp_cmp )

	for eit in eitems:
		item_cur = eit[0]
		for item in eit[1:]:
			parents[item[0]] = { item_cur[0] : (0,"e") }

	remaining = filterdict( remaining , parents )

	print "done filtering equal vectors"

	for m in range(1,3):
		if len( remaining ) <= m:
			break
		print "filtering linear dependence of order " , m
		Ls = rp_line_finder( remaining , m )
		gg = gen_graph( remaining , Ls , m )
		for u in gg:
			if gg[u]:
				parents[u] = dict( [ (v,(m,"lc")) for v in gg[u] ] )

		remaining = filterdict( remaining , parents )

		print "%s vectors remaining" % ( len(remaining,) )

		ggdot = graph.Dot( gg )
		ggdot.save_image( "foo%d.png" % (m,) , "png")

	for r in remaining:
		parents[r] = {}

	return parents


def main_gg():
	shape="triangle"
	degree=3
	A0dict=build_tensors.laplacian(shape,degree)
	process_gen_graph(A0dict)

def main():
    import string
    results = {}
    shape = "triangle"
    degrees = range(2,4)
    max_level = 4
    for degree in degrees:
        A0dict = \
               build_tensors.weighted_laplacian_coeff_first(shape,degree)
	results[degree] = process2(A0dict,max_level)

    # now let's print the table
    # first row is heading
    result_str = """
\\begin{tabular}{%s}
%s \\\\ \\hline
%s \\\\
%s \\\\
%s \\\\
%s \\\\
%s \\\\ \\hline
""" % ( string.join( ["c"]*(len(degrees)+1) , "|" ) , \
	string.join( [""] + map(str,degrees) , " & ") , \
	string.join( ["$n$"] + [str(a[1]) \
				for a in results.values()] , " & ") , \
	string.join( ["$d$"] + [str(a[2]) \
				for a in results.values()] , " & ") , \
	string.join( ["num zero"] + [str(a[3]) \
			       for a in results.values()] , " & ") , \
	string.join( ["num equal"] + [str(a[4]) \
			       for a in results.values()] , " & ") , \
	string.join( ["num colinear"] + [str(a[6][0]) \
			       for a in results.values()] , " & " ) )

    for level in range(2,max_level+1):
	result_str += """%s \\\\
%s \\\\
%s \\\\
%s \\\\ \\hline
""" % ( string.join( ["size for %d-search"%(level,)] + \
		     [ str(a[5][level-1]) \
		       for a in results.values()] , " & " ) , \
	string.join( ["num with %d-dependency" % (level,) ] + \
		     [ str(a[6][level-1]) \
		       for a in results.values()] , " & " ) , \
	string.join( ["generator size"] + \
		     [ str(a[5][level]) \
		       for a in results.values()] , " & " ) , \
	string.join( ["MAPs"] + \
		     [ str(a[7][level-1]) \
		       for a in results.values()] , " & " ) )


    result_str += "\\end{tabular}"

    print result_str

if __name__ == "__main__":
    main_gg()