Vanishing component analysis code in Python, using sympy

result.txt

[-0.999999999990622*x0 + 3.93039236601389e-7*x1 + 4.31294626702504e-6*x2 + 0.000105056139124062, 4.31294796214237e-6*x0 + 4.31294626702504e-6*x1 + 0.999999999981398*x2 - 1.00012486669319, 7.72281068263184e-14*x0**2 + 3.92997720318555e-7*x0*x1 - 1.6949787169904e-12*x0*x2 - 4.30611426233118e-11*x0 + 0.499970847799378*x1**2 - 4.31269649833429e-6*x1*x2 - 0.000109564584584033*x1 + 9.30021778700191e-12*x2**2 + 4.7254635180911e-10*x2 - 0.499999993997451]
[-0.707104306393211*x0 + 1.46301599452756e-5*x1 + 0.707109255819873*x2 - 0.707064928651531, 0.707109255878223*x0 - 1.58849047682263e-6*x1 + 0.707104306484427*x2 - 0.707161709963975, 6.57532422227314e-11*x0**2 + 1.14669699610752e-5*x0*x1 - 1.05747170661535e-10*x0*x2 - 8.88539783667339e-10*x0 + 0.499942647857534*x1**2 - 9.22083526571183e-6*x1*x2 - 7.747799700903e-5*x1 + 4.25167783553563e-11*x2**2 + 7.14493802638299e-10*x2 - 0.499999996998236]
[-0.164419808895816*x0 - 2.16240662386745e-5*x1 + 0.986390453104177*x2 - 0.986462928050236, 0.986390453341201*x0 - 3.67266170431257e-6*x1 + 0.164419808854811*x2 - 0.164520223495808, 2.26170212244376e-15*x0**2 + 6.72589109738161e-8*x0*x1 + 1.47523212407374e-12*x0*x2 - 4.47708740776069e-12*x0 + 0.500039446009784*x1**2 + 2.19353612295282e-5*x1*x2 - 6.65702216910878e-5*x1 + 2.40561057787267e-10*x2**2 - 1.46012666758216e-9*x2 - 0.499999997784377]

VCAsample.py

# -*- coding: utf-8 -*-
"""
Created on Sat Apr 12 21:02:25 2014

@author: xiangze
"""

import numpy as np
from sympy import *
from scipy import linalg, mat, dot

def limits(f,l):
    v=f
    for i in l:
        v=limit(v,i[0],i[1])
    return v

def limitsv(f,var,val):
    return [limits(f,zip(var,v)) for v in val]

"""
FindRangeNull (using SVD)
input
F: non-vanishing polynomials
C: Canditate polynomials
S: zero points 
vars: list of variables 
eps: minimum eigen value

return
V: vanishing polynomials
F: non-vanishing polynomials
"""

def FindRangeNull(F,C,S,vars,eps):
    k=len(C)
    At=[]
    Fs=[]

    for i in range(k):
        fi=limitsv(C[i],vars,S)
        f=C[i]-sum([np.dot(fi,limitsv(g,vars,S))*g for g in F])
        Fs.append(f)
        At.append(limitsv(f,vars,S))

    """
    print "S",S
    print "vars",vars
    print "F",F
    print "Fs",Fs
    print "At",mat(At)
    """

    L, D, U = linalg.svd(mat(At).transpose())

    """
    print "num:k",k
    print "At",mat(At).shape, "L,D,U", L.shape, D.shape, U.shape
    print "D:",D
    """

    G=[sum([U[i,j]*Fs[j]for j in range(len(Fs))]) for i in range(k)]
    G=[expand(i) for i in G]
#    print "G:",G

    V=list()
    F=list()
    for i in range(k):
        if( i<len(D) and D[i]>eps):
            Fv=limitsv(G[i],vars,S)
#            print "Fvalue",Fv
#            p=np.linalg.norm(Fv)
            p=sum([n*n for n in Fv])**0.5
            F.append(G[i]/p)
        else:
            V.append(G[i])
    return F,V

def VCAsample(S,eps=0.01,maxite=100):
    n=len(S[0])
    m=len(S)
    vars = [Symbol('x'+str(i))for i in range(n)]
    F=[1/(m**0.5)]
    C=vars

#    print "init:",F

    (F1,V1)=FindRangeNull(F,C,S,vars,eps)
    F=F+F1
    V=V1
    """
    print "C",C
    print "F",F
    print "V",V
    """

    Ft=F1
    ind=1
    while(ind<maxite):
#        print "iteration:",ind
        C=[simplify(g*h) for g in Ft for h in F1 ]
        if not C:
            break
#        print "C",len(C),";",C
        (Ft,Vt)=FindRangeNull(F,C,S,vars,eps)
        F=F+Ft
        V=V+Vt
        ind=ind+1

    return F,V

def test_randVCA(dataorg,eps=0.002,n=3):
    for seed in range(n):
        np.random.seed(seed)
        data=[ np.random.rand(3)*eps+i for i in dataorg]
        (FA,VA)=VCAsample(data,eps,400)
        print VA
        #return  VA
    
def test_VCAsample():
    N1=3.
    N2=3.
    Nth=2.
    eps=0.0002

    zs=[i/N1 for i in range(0,int(N1))]
    z1s=[i/N2 for i in range(0,int(N2))]
    ths=[2*pi*i/Nth for i in range(0,int(Nth))]

    dataorg=[(sin(th),cos(th),1) for th in ths ]
    test_randVCA(dataorg,eps)

    #dataAorg=[(sin(th),cos(th),z) for th in ths for z in zs ]
    #test_randVCA(dataAorg,eps)
    
    #dataBorg=[((1-z**2)*sin(th),(1-z**2)*100*cos(th),z) for th in ths for z in z1s ]
    #V=test_randVCA(dataBorg,eps)    


if __name__ == "__main__":
    test_VCAsample()