helpers.py

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

# -*- coding: iso-8859-1 -*-

#
# This file is a part of the Bayes Blocks library
#
# Copyright (C) 2001-2006 Markus Harva, Antti Honkela, Alexander
# Ilin, Tapani Raiko, Harri Valpola and Tomas 謘tman.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2, or (at your option)
# any later version.
#
# This program 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 General Public License (included in file License.txt in the
# program package) for more details.
#
# $Id: Helpers.py 7 2006-10-26 10:26:41Z ah $
#

from Label import *
import Net
try:
    import numpy.oldnumeric as Numeric
    import numpy as MLab
    import numpy.random as RandomArray
    from numpy.linalg import inv, eigh
except:
    import Numeric
    from LinearAlgebra import inverse as inv
    from LinearAlgebra import Heigenvectors as eigh
    import RandomArray
    import MLab
import string
import fileinput
import re

def LoadMatlabArray(file, array):
    try:
        mat = Net.MatArray()
    except:
        raise RuntimeError("Matlab arrays are not supported in this library version")
    mat.LoadFromFile(file, array)
    res = []
    for i in range(mat.Size(2)):
        res.append(mat.ExtractDV(i))
    return res

def CostDifference(node):
    grad = Net.DSSet()
    val = Net.DSSet()
    node.ChildGradReal(grad)
    if (grad.ex):
        node.GetReal(val, Net.DFlags(1,1,1))
    else:
        node.GetReal(val, Net.DFlags(1,1))
        val.ex = 1
    cdiff = node.Cost() + grad.var * (val.var - val.mean * val.mean) +\
            grad.mean * val.mean + grad.ex * (val.ex - 1)
    return -cdiff

def CostDifferenceV(node):
    grad = Net.DVSet()
    val = Net.DVH()
    node.ChildGradRealV(grad)
    if (len(grad.ex)):
        node.GetRealV(val, Net.DFlags(1,1,1))
    else:
        node.GetRealV(val, Net.DFlags(1,1))
    cdiff = node.Cost()
    for i in range(len(grad.mean)):
        cdiff += grad.mean[i] * val.Mean(i)
    for i in range(len(grad.var)):
        cdiff += grad.var[i] * (val.Var(i) - val.Mean(i) * val.Mean(i))
    for i in range(len(grad.ex)):
      cdiff += grad.ex[i] * (val.Exp(i) - 1)
    return -cdiff

def LoadAsciiData(file):
    res = []
    for line in fileinput.input(file):
        s = map(float, string.split(line))
        tmp = Net.DV(len(s))
        for i in range(len(s)):
            tmp[i] = s[i]
        res.append(tmp)
    return res

def ParentList(node):
    res = []
    i = 0
    newpar = node.GetParent(i)
    while newpar:
        res.append(newpar)
        i += 1
        newpar = node.GetParent(i)
    return res

def VariableParents(node):
    nodes=[node]
    variables=[]
    while(len(nodes)):
        for n in ParentList(nodes.pop()):
            if isinstance(n, PyNet.Variable):
                variables.append(n)
            else:
                nodes.append(n)
    return variables

def ChildList(node):
    res = []
    i = 0
    newpar = node.GetChild(i)
    while newpar:
        res.append(newpar)
        i += 1
        newpar = node.GetChild(i)
    return res

def FindLeaf(rootnode, parind):
    parent = rootnode.GetParent(parind)
    if not re.match('Sum2', parent.GetType()):
        return [rootnode, parent]
    roots = [parent]
    while 1:
        newroots = []
        for root in roots:
            for par in ParentList(root):
                if re.match('Sum2', par.GetType()):
                    newroots.append(par)
                else:
                    return [root, par]
        roots = newroots

def AddSum2(net, rootnode, parind, newpar, label):
    leaf = FindLeaf(rootnode, parind)
    tmp = Net.Sum2(net, label, leaf[1], newpar)
    leaf[0].ReplacePtr(leaf[1], tmp)
    tmp.AddChild(leaf[0])
    leaf[1].NotifyDeath(leaf[0])

def AddSum2V(net, rootnode, parind, newpar, label):
    leaf = FindLeaf(rootnode, parind)
    tmp = Net.Sum2V(net, label, leaf[1], newpar)
    leaf[0].ReplacePtr(leaf[1], tmp)
    tmp.AddChild(leaf[0])
    leaf[1].NotifyDeath(leaf[0])

def AddWeight(net, inode, onode, weightlabel, mnode, vnode,
              addcost = 0, labelbase = ""):
    weight = Net.Gaussian(net, weightlabel, mnode, vnode)
    tmp = Net.ProdV(net, Label(labelbase + "_newprod"), weight, inode)
    AddSum2V(net, onode, 0, tmp, Label(labelbase + "_newsum"))
    weight.Update()
    if not net.TryPruning(weight, 0, addcost):
        return weight

def GetMean(node):
    if IsSequence(node):
        return Numeric.array([GetMean(n) for n in node], Numeric.Float)
    else:
        if node is None:
            return 0.0
        else:
            tmp = Net.DSSet()
            node.GetReal(tmp, Net.DFlags(1))
            return tmp.mean

def GetVar(node):
    if IsSequence(node):
        return Numeric.array([GetVar(n) for n in node], Numeric.Float)
    else:
        if node is None:
            return -1
        else:
            tmp = Net.DSSet()
            node.GetReal(tmp, Net.DFlags(0, 1))
            return tmp.var

def GetExp(node):
    if IsSequence(node):
        return Numeric.array([GetExp(n) for n in node], Numeric.Float)
    else:
        if node is None:
            return -1
        else:
            tmp = Net.DSSet()
            node.GetReal(tmp, Net.DFlags(0, 0, 1))
            return tmp.ex

def GetMeanV(node):
    if IsSequence(node):
        return Numeric.array([GetMeanV(n) for n in node], Numeric.Float)
    else:
        tmp = Net.DVH()
        if not node.GetRealV(tmp, Net.DFlags(1)):
            return
        if tmp.vec:
            return Numeric.array(tmp.vec.mean, Numeric.Float)
        else:
            return tmp.scalar.mean

def GetVarV(node):
    if IsSequence(node):
        return Numeric.array([GetVarV(n) for n in node], Numeric.Float)
    else:
        tmp = Net.DVH()
        if not node.GetRealV(tmp, Net.DFlags(0, 1)):
            return
        if tmp.vec:
            return Numeric.array(tmp.vec.var, Numeric.Float)
        else:
            return tmp.scalar.var

def GetExpV(node):
    if IsSequence(node):
        return Numeric.array([GetExpV(n) for n in node], Numeric.Float)
    else:
        tmp = Net.DVH()
        if not node.GetRealV(tmp, Net.DFlags(0, 0, 1)):
            return
        if tmp.vec:
            return Numeric.array(tmp.vec.ex, Numeric.Float)
        else:
            return tmp.scalar.ex

def GetDiscrete(node):
    if IsSequence(node):
        return Numeric.array([GetDiscrete(n) for n in node], Numeric.Float)
    else:
        return Numeric.array(node.GetDiscrete(), Numeric.Float)

def GetDiscreteV(node):
    if IsSequence(node):
        return Numeric.array([GetDiscreteV(n) for n in node], Numeric.Float)
    else:
        tmp = Net.VDDH()
        node.GetDiscreteV(tmp)
        tdim = tmp.vec.size()
        kdim = tmp.vec.DDsize()
        ary = Numeric.zeros((kdim, tdim), Numeric.Float)
        for i in range(kdim):
            for j in range(tdim):
                val = tmp.vec.Get(j, i)
                ary[i,j] = val
        return ary

def GetLabel(node):
    if IsSequence(node):
        return map(GetLabel, node)
    else:
        return node.GetLabel()

def GetType(node):
    if IsSequence(node):
        return map(GetType, node)
    else:
        return node.GetType()

def GetCost(variable):
    if type(variable) == type([]):
        return map(GetCost, variable)
    else:
        return variable.Cost()

def DoPCA(x, dim, returnA=0, returnW=0, returnd=0, removemean=1):
    xdim = x.shape
    if removemean:
        x0 = x - MLab.mean(x, 1)[:, Numeric.NewAxis]