helpers.py

The library is a C++/Python implementation of the variational building block framework introduced in

    else:
        x0 = x
    C = Numeric.dot(x0, Numeric.transpose(x0)) / xdim[1]
    d, V = eigh(C)
    I = Numeric.argsort(-d)
    I2 = list(I[0:dim])
    V2 = Numeric.take(V, I2)
    d2 = Numeric.take(d, I2) ** -.5
    W = Numeric.dot(MLab.diag(d2), V2)
    y = Numeric.dot(W, x0)
    val = [y]
    if returnA:
        A = Numeric.dot(Numeric.transpose(V2),
                        MLab.diag(Numeric.take(d, I2) ** .5))
        val.append(A)
    if returnW:
        val.append(W)
    if returnd:
        val.append(Numeric.take(d, I2))
    if len(val) == 1:
        return val[0]
    else:
        return tuple(val)

def NonNegPCA(x, dim, iters=1000, orthiter=10,
              epsilon=1e-10, verbose=0, printiter=50, returnW=0,
              learningscheme={0: 1e-4, 10: 1e-3}):
    """Non-Negative PCA. Suited for positive sources. Despite its
    name performs positive ICA."""
    def mnorm(x):
        return Numeric.sqrt(Numeric.sum(Numeric.sum(x**2))) / (
            x.shape[0] * x.shape[1])
    trans = Numeric.transpose
    mmul = Numeric.dot
    
    z = DoPCA(x, dim, removemean=0)
    z = z / MLab.std(z, 1)[:,Numeric.NewAxis]
    W = Numeric.identity(z.shape[0], Numeric.Float)
    dW = Numeric.array(W, Numeric.Float)
    iter = 0
    rate = learningscheme[0]

    while iter < iters and mnorm(dW) > epsilon:
        y = mmul(W, z)
        fy = y * (y < 0)
        dW = -mmul(mmul(fy, trans(y)) - mmul(y, trans(fy)), W)
        W += rate * dW
        if iter % orthiter == 0:
            W = trans(Orthogonalize(trans(W)))
        if verbose and iter % printiter == 0:
            print "%d : ||dW|| = %f" % (iter, mnorm(dW))
        if iter in learningscheme:
            rate = learningscheme[iter]
        iter += 1

    y = mmul(W, z)
    if returnW:
        return y, W
    else:
        return y

def Orthogonalize(basis):
    """Orthogonalizes the set of column vectors of a matrix."""
    def norm(v):
        return Numeric.sqrt(Numeric.sum(v**2))

    q = Numeric.zeros(basis.shape, Numeric.Float)
    q[:,0] = basis[:,0] / norm(basis[:,0])

    for j in range(1, basis.shape[1]):
        vcol = basis[:,j:j+1]
        vrow = basis[:,j]
        u = vrow - Numeric.sum(Numeric.sum(vcol * q[:,:j]) * q[:,:j], -1)
        q[:,j] = u / norm(u)

    return q


def EmbeddedPCA(data, sdim, emblen, returnA = 0, returnB=0):
    xdim, tdim = data.shape
    if type(emblen) != type(()):
        emblen = (emblen, emblen)
    edata = Numeric.zeros(((1+emblen[0]+emblen[1])*xdim,
                           tdim + (emblen[0]+emblen[1]))).astype('d')
    for i in range(-emblen[0], emblen[1]+1):
        edata[(i+emblen[0])*xdim:(i+emblen[0]+1)*xdim,
              (i+emblen[0]):(tdim+i+emblen[0])] = data[:, :]

    C = MLab.cov(edata[:, (emblen[0]+emblen[1]):tdim], rowvar=1)
    d, V = eigh(C)
    I = Numeric.argsort(-d)
    I2 = list(I[0:sdim])
    V2 = Numeric.take(V, I2)
    d2 = Numeric.take(d, I2) ** -.5
    d3 = Numeric.take(d, I2) ** .5
    y = Numeric.dot(Numeric.dot(MLab.diag(d2), V2),
                    (edata[:, emblen[0]:(tdim+emblen[0])]
                     - MLab.mean(edata, 1)[:, Numeric.NewAxis]))
    B = Numeric.dot(Numeric.transpose(V2), MLab.diag(d3))

    A = Numeric.dot(
        Numeric.dot(y[:, 1:], Numeric.transpose(y[:, :-1])),
        inv(Numeric.dot(
        y[:, :-1], Numeric.transpose(y[:, :-1]))))

    if returnB:
        if returnA:
            return y, A, B
        else:
            return y, B
    else:
        return y

def SBSS(X, dim, f, iters = 1000, epsilon = 1e-4,
         verbose = 0, printiter = 50):
    """Semi-blind source separation. f is some noisereduction function
    which should take a one-dimensional array as an argument and return
    one too."""
    def norm(x):
        return Numeric.sqrt(Numeric.sum(x**2))

    mmul = Numeric.matrixmultiply
    trans = Numeric.transpose

    xdim = X.shape[0]
    tdim = X.shape[1]

    X = DoPCA(X, xdim)
    w = RandomArray.standard_normal(xdim)
    W = Numeric.zeros((dim, xdim), Numeric.Float)
    S = Numeric.zeros((dim, tdim), Numeric.Float)
    Xplus = trans(mmul(inv(mmul(X, trans(X))), X))

    for i in range(dim):
        if verbose:
            print "Estimating signal %d" % i
        wold = Numeric.zeros(xdim)
        iter = 0
        while norm(wold - w) > epsilon and iter < iters:
            wold = w
            s = mmul(w, X)
            s = f(s)
            w = mmul(s, Xplus)
            for j in range(i):
                w = w - mmul(W[j,:], w) * W[j,:]
            w = w / norm(w)
            iter += 1
            if verbose and iter % printiter == 0:
                print "%d : ||wold - w|| = %f" % (iter, norm(wold - w))
        W[i,:] = w
        S[i,:] = mmul(w, X)

    return S

def Array2DV(ary):
    ary = Numeric.array(ary) # now it is certainly an array
    #assert len(ary.shape) == 1
    # with the typemaps defined, we can just convert the array to a list
    # and everything should work just fine
    return ary.tolist()

def Array2VDD(a):
    if len(a.shape) != 2:
        raise ValueError, "Expected a two dimensional array"
    v = Net.VDD(a.shape[1], a.shape[0])
    for i in range(a.shape[1]):
        for j in range(a.shape[0]):
            v.Set(i, j, a[j,i])
    return v

def Array2IntV(ary):
    ary = Numeric.array(ary) # now it is certainly an array
    #assert len(ary.shape) == 1
    # with the typemaps defined, we can just convert the array to a list
    # and everything should work just fine
    return ary.tolist()

def GetDescendants(node):
    checked = {}
    if type(node) in (type(()), type([])):
        check = {}
        for i in node:
            check[i.GetLabel()] = i
    else:
        check = {node.GetLabel():node}
    while len(check) > 0:
        k = check.popitem()[1]
        for p in ChildList(k):
            if not checked.has_key(p.GetLabel()):
                check[p.GetLabel()] = p
                checked[p.GetLabel()] = p
    return checked.values()

def GetVariableDescendants(net, node):
    descendants = GetDescendants(node)
    variable_nodes = []
    for i in descendants:
        tmp = net.GetVariable(i.GetLabel())
        if tmp is not None:
            variable_nodes.append(tmp)
    return variable_nodes

def ApplyTestData(net, inputnode, inputrange, numiter = 5,
                  printcost = 0, variables = []):
    tdim = net.Time()
    if type(inputnode) not in (type(()), type([])):
        inputnode = (inputnode,)
    if len(inputnode) not in (1,2):
        raise ValueError, "Length of inputnode must be 1 or 2"
    if type(inputrange[0]) not in (type(()), type([])):
        inputrange = (inputrange,)
    if len(inputrange) != len(inputnode):
        raise ValueError, "Length of inputrange must equal to lengt of inputnode"
    for i in range(len(inputrange)):
        if len(inputrange[i]) != 2:
            raise ValueError, "Length of inputranges elements must be 2"
    
    for v in variables:
        v.Unclamp()

    if len(inputnode) == 2:
        r = int(math.sqrt(tdim))
        x = Numeric.fromfunction(lambda x,y:x, (r,r))
        y = Numeric.fromfunction(lambda x,y:y, (r,r))
        x = x*(inputrange[0][1]-inputrange[0][0])/(r-1)+inputrange[0][0]
        y = y*(inputrange[1][1]-inputrange[1][0])/(r-1)+inputrange[1][0]
        x = Numeric.resize(x, (tdim,))
        y = Numeric.resize(y, (tdim,))
        inputnode[0].Clamp(Array2DV(x))
        inputnode[1].Clamp(Array2DV(y))
    else:
        #Clamps lienarly spaced values in inputrange to node
        inputnode[0].Clamp(Array2DV(
            Numeric.arrayrange(tdim)*(inputrange[0][1]-inputrange[0][0]) /
            (tdim-1)+inputrange[0][0]))

    iter = 0
    while iter < numiter:
        map(Net.Variable.Update, variables)
        iter += 1
        if printcost and (iter%printcost == 0):
            print iter, ":", net.Cost()
    if len(inputnode) == 2:
        return (r,r)
    else:
        return (tdim,)

def IsInf(val):
    return (val > 1e10) and (val/2 == val)

def IsNaN(val):
    return ((val > val) or ((0 == val) and (1 == val)))

def IsSequence(val):
    return (type(val) in [type([]), type(()), Numeric.ArrayType])

def __UniqueHelper(x, y):
    if x[-1]==y:
        return x
    else:
        return x+[y]

def Unique(l):
    """Return a list of unique elements in a list, i.e. leave only one of
    equal consecutive elements."""
    if len(l) < 2:
        return l
    else:
        return reduce(__UniqueHelper, l[1:], [l[0]])